This paper focuses on fault-related uncertainties in the subsurface, which can significantly affect the numerical simulation of physical processes. Our goal is to use dynamic data and process-based simulation to update structural uncertainty in a Bayesian inverse approach. We propose a stochastic fault model where the number and features of faults are made variable. In particular, this model samples uncertainties about connectivity between the faults. The stochastic three dimensional fault model is integrated within a stochastic inversion scheme in order to reduce uncertainties about fault characteristics and fault zone layout, by minimizing the mismatch between observed and simulated data. The stochastic fault model uses a priori information such as fault orientation, location, size and sinuosity, to sample both geometrical and topological uncertainties with realistic fault descriptions. Each fault object is parameterized by the random vector used to simulate fault features. Then, during inversion, the random vector of the current model is stochastically perturbed, producing a new parameter vector used as input by the stochastic fault model to produce a new model. Even if the topology varies from one model to another, the algorithm produces correlated models so that their flow responses evolve quite smoothly. The methodology is applicable in general and illustrated on a synthetic two-phase flow example. A first set of models is generated to sample the prior uncertainty space. Then models minimizing reference water-saturation data misfit are used as seeds to generate continuous Monte Carlo Markov Chains (MCMC) of models with discrete states. Posterior models reduce uncertainties about fault position, while the topology varies from one model to another. A second example highlights the interest of the parameterization when interpreted data are available, by perturbing geological scenarios and falsifying those that do not match two-phase flow observations.