This paper presents a novel methodology for inverse modeling of groundwater flow and transport problems in a Monte Carlo framework, i.e., multiple solutions to the inverse problem are generated. The methodology is based on the concept of random mixing of spatial random fields. The conditional target hydraulic transmissivity field is obtained as a linear combination of unconditional spatial random fields. The corresponding weights of the linear combination are selected such that the spatial variability of the hydraulic transmissivities as well as the actual observed transmissivity values are reproduced. The constraints related to the hydraulic head and contaminant concentration observations are nonlinear. In order to fulfill these constraints, a specific property of the presented approach is used. A connected domain of fields fulfilling all linear constraints is identified. This domain includes an infinite number of realizations, and in this domain, the head and concentration deviations are minimized using standard continuous optimization techniques. The methodology uses spatial copulas to describe the spatial dependence structure. A combination with multiple point statistics allows inversion under specific structural constraints.