History matching is an integral part of field development planning for oil and gas reservoirs. It can be viewed as an inverse modeling technique which utilizes information contained in the observed flow response variables like flow rate, well bottom hole pressure etc. to better quantify the spatial distribution of reservoir model parameters like permeability or porosity. As the reservoir parameters and flow response variables are often related by non-linear relationships, the solution to the inverse problem of history matching is often non-unique. This makes it a problem that can be better handled by stochastic approaches than deterministic approaches. The proposed approach called 'Indicator-based Data Assimilation' (InDA) is suitable for such problems. To initiate the process, multiple realizations of the reservoir model parameters are generated based on the available information related to the reservoir description. The ensemble of realizations serves as a measure of initial uncertainty in the spatial distribution of model parameters. Next, flow simulations are run on the ensemble of models and the difference between the observed and simulated flow response variables are used to update the reservoir parameters using InDA. With successive updates, the uncertainty in the model parameters is reduced and the spatial distribution approaches the "true" distribution. Considering the residual uncertainty of the final updated models, reliable field development planning decisions can be made. The proposed method is validated using a realistic reservoir with complex channel-like features emulating a reservoir formed in a fluvial depositional environment. Liquid rate data from existing wells are used for updating reservoir parameters for several time steps using InDA. A comparison of the spatial distribution of final model parameters with the reference model used for validation shows a good match. After the history matching period, existing and new infill wells are run in a forecast mode where the observed and simulated flow responses show a good match. As majority of oil reservoirs comprise of high permeability oil-bearing zones in form of channels passing across low permeability zones, the statistical permeability distribution is bimodal making it non-Gaussian. It is shown that 'Indicator Transformation' of variables used in InDA preserves the non-Gaussian structure of the permeability field in comparison to methods like 'Ensemble Kalman Filter' (EnKF) that are sub-optimal in such cases.