Accurately characterizing fractures is complex. Several studies have proposed reducing uncertainty by incorporating fracture characterization into simulations, using a probabilistic approach, to maintain the geological consistency, of a range of models instead of a single matched model. We propose a new methodology, based on one of the steps of a general history-matching workflow, to reduce uncertainty of reservoir attributes in naturally fractured reservoirs. This methodology maintains geological consistency and can treat many reservoir attributes. To guarantee geological consistency, the geostatistical attributes (e.g., fracture aperture, length, and orientation) are used as parameters in the history matching. This allows us to control Discrete Fracture Network attributes, and systematically modify fractures. The iterative sensitivity analysis allows the inclusion of many (30 or more) uncertain attributes that might occur in a practical case. At each uncertainty reduction step, we use a sensitivity analysis to identify the most influential attributes to treat in each step. Working from the general history-matching workflow of Avansi et al. (2016), we adapted steps for use with our methodology, integrating the history matching with geostatistical modeling of fractures and other properties in a big loop approach. We applied our methodology to a synthetic case study of a naturally fractured reservoir, based on a real semi-synthetic carbonate field, offshore Brazil, to demonstrate the applicability in practical and complex cases. From the initial 18 uncertain attributes, we worked with only 5 and reduced the overall variability of the Objective Functions. Although the focus is on naturally fractured reservoirs, the proposed methodology can be applied to any type of reservoir.