Modeling pixel intensity with the normal distribution
The first thing to consider is that, even continuous measurements are actually discrete. For instance, temperature is always considered as continuous but every time you process temperature data, these are discrete data to the precision of the sensor used to measure them. Same is true for every physical, measurable quantity (including intensity). That being said, with our modern physics models, it makes no sense to consider that temperatures follow discrete mass probabilities with several thousands possible output and everyone studies them under continuous assumptions. 

So where do we set the boundary? When do we consider that there are "too many" possible outputs to work with probability mass functions and choose to work with continuous distributions? This is an open question! :)

Here is the thing, we do not lose information when processing discrete data in a continuous way so there is no restriction for doing so!

Now consider the pixel intensities of an image, most of the time, we would need to consider discrete mass probabilities with 256 possible outputs, meaning 256 parameters (or weights) to be estimated or tuned. If we model the intensity by a continuous distribution, the number of parameters is drastically reduced! For a Gaussian for instance: 2 parameters. Even using more complex distributions or mixture models, we will not need to tune hundreds of parameters. This is the main advantage! Compact representation!

About the choice of the Gaussian specifically, let say that it has been empirically proven to work great. The Central Limit Theorem also supports this common choice. And now, on the top of this, many algorithms are readily available under the assumption that the data follows Gaussian distributions making them even more used by people in the field. Doesn't mean that nothing else could be used but this is the reality of things.