Abstract Fractured reservoirs are highly heterogeneous and can be characterized by the probability distributions of fracture properties in a discrete fracture network model. The relationship between production performance and the fracture parameters is vastly nonlinear, rendering the process of adjusting model parameters to match both the static geological and dynamic production data challenging. This creates a need for a comprehensive history matching workflow for fractured reservoirs, which considers different local as well as global fracture parameters and leads to multiple equally-probable realizations of the discrete fracture network model parameters for uncertainty quantification. This paper presents an integrated approach for the history matching of fractured reservoirs. This new methodology includes generating multiple discrete fracture models, upscaling them for numerical multiphase flow simulation, and updating the fracture properties using dynamic flow responses such as continuous rate and pressure measurements. Available geological and tectonic information such as well-logs, seismic, and structural maps are incorporated into commercially available DFN modeling and simulation software to infer the probability distributions of relevant fracture parameters (including aperture, length, connectivity, and intensity) and to generate multiple discrete fracture network model realizations. The fracture models are further upscaled into an equivalent continuum dual-porosity model in the software using either analytical approaches or dynamic methods. The upscaled models are subjected to the flow simulation, and their production performances are compared to the true recorded responses. An automated history matching algorithm is implemented to reduce the uncertainties of the fracture properties. Components of vectors representing the principal flow directions and average fracture orientations are obtained by means of eigenvector decomposition of the permeability tensor and are optimized in the algorithm. In addition, both global fracture intensity and the local grid based intensity, which highly affect the fluid flow pattern and rate in different regions of the reservoir, are adjusted. A case study with various fracture sets is presented. The initial realizations were generated by means of Monte Carlo simulations, using the observed fractures at the well locations. Fracture intensity, orientation, and conductivity of different fracture sets were the uncertain parameters in our studies. Using the proposed methodology, parameters of different fracture sets were satisfactorily updated. Implementation of this automated history matching approach resulted in multiple equally probable discrete fracture network models and their upscaled flow simulation models that honor the geological information, and at the same time they match the dynamic production history.