Inverse modeling aims to adjust the uncertain parameters of the geological model so that the numerical simulation results can reproduce the observed data and then provide reliable predictions, which is of great significance for the development of subsurface resources. Nowadays, inverse modeling is no longer limited to obtaining a single posterior model but evaluates the uncertainty of the predictions in terms of the bandwidth of dynamic responses by applying Monte Carlo methods to obtain a suite of equiprobable geological models (Arnold et al., 2013;Suzuki & Caers, 2006). Among the ensemble-based methods, the ensemble Kalman filter (EnKF; Aanonsen et al., 2009;Agbalaka & Oliver, 2008) is one of the most successful and effective methods. However, the continuous restart of the simulator may prevent the use of EnKF (Emerick & Reynolds, 2013a). Unlike EnKF, the ensemble smoother (ES; Van Leeuwen & Evensen, 1996) can avoid this problem by calculating a global update with all available observations. Nevertheless, the single update also hinders adequate data matching. More effective methods are proposed to overcome this problem by implementing the iterative forms of ES, such as the ES with multiple data assimilation (ESMDA) (Emerick & Reynolds, 2013a, 2013b). However, one intrinsic limitation is that these data assimilation methods are based on Gaussian assumptions, which makes them more suitable for geological parameters with the continuous spatial distribution. The posterior geological model parameters obtained by the data