Simulation of random fields is widely used in Earth sciences for modeling and uncertainty quantification. The spatial features of these fields may have a strong impact on the forecasts made using these fields. For instance, in flow and transport problems the connectivity of the permeability fields is a crucial aspect. Multi-Gaussian random fields are the most common tools to analyze and model continuous fields. Their spatial correlation structure is described by a covariance or variogram model. However, these types of spatial models are unable to represent highly or poorly connected structures even if a broad range of covariance models can be employed. With this type of model, the regions with values close to the mean are always well connected whereas the regions of low or high values are isolated. Substitution random functions (SRFs) belong to another broad class of random functions that are more flexible. SRFs are constructed by composing (\(Z=Y\circ T\)) two stochastic processes: the directing function T (latent field) and the coding process Y (modifying the latent field in a stochastic manner). In this paper, we study the properties of SRFs obtained by combining stationary multi-Gaussian random fields for both T and Y with bounded variograms. The resulting SRFs Z are stationary, but as T has a finite variance Z is not ergodic for the mean and the covariance. This means that single realizations behave differently from each other. We propose a simple technique to control which values (low, intermediate, or high) are connected. It consists of adding a control point on the process Y to guide every single realization. The conditioning to local values is obtained using a Gibbs sampler.