Multiple-point statistical (MPS) simulation methods have gained widespread adoption across various Earth science disciplines. They offer a versatile framework for simulating intricate spatial patterns and heterogeneity in both surface and subsurface structures. While these simulations adeptly incorporate conditioning to hard data, such as information from boreholes, conditioning to indirect data (e.g. geophysical data) is more challenging. A new methodology is introduced that provides geostatistical realisations honouring indirect geophysical data and complex prior knowledge described by a training image. An MPS simulation is iteratively built up pixel-by-pixel starting from an empty grid or with initial hard conditioning data if available. During each simulation step, a pixel value is selected from a set of candidates proposed by the MPS algorithm. This selection is made proportionally to an approximated likelihood that accounts for indirect geophysical data. The expected values and uncertainty quantification are obtained by simulating many complete field realisations. Our approach, which we name Indirect Data Conditional Simulations (IDCS), is tested for multi-Gaussian and complex subsurface structures with synthetic data from linear and non-linear crosshole ground-penetrating radar responses. The IDCS method is inherently approximate due to the finiteness of the training image, a limited number of MPS candidates at each simulation step and the need to approximate intractable likelihood functions. Nevertheless, the results demonstrate that the posterior approximations obtained by IDCS are often comparable to those obtained with a Markov chain Monte Carlo method, with IDCS being at least one order of magnitude faster. While the method performs the best when the underlying physics is modelled as a linear response, encouraging preliminary results considering non-linear physical responses are provided.