This paper studies a Gaussian process approximation for a class of stochastic traffic flow models. It can be used to efficiently and accurately evaluate the joint (in the spatial and temporal sense) distribution of vehicle-density distributions in road traffic networks of arbitrary topology. The Gaussian approximation follows, via a scaling-limit argument, from a Markovian model that is consistent with discrete-space kinematic wave models. We describe in detail how this formal result can be converted into a computational procedure. The performance of our approach is demonstrated through a series of experiments that feature various realistic scenarios. Moreover, we discuss the computational complexity of our approach by assessing how computation times depend on the network size. We also argue that the (debatable) assumption that the vehicles’ headways are exponentially distributed does not negatively impact the accuracy of our approximation.