Simulating fluid flow in complex fracture systems encountered in subsurface systems, such as those encountered in tight or shale gas reservoirs, presents considerable challenges. Sophisticated numerical techniques and characterization approaches are required to integrate static and dynamic data from multiple sources to account for the intricate interplay between fractures and the surrounding rock matrix. A novel workflow using an indicator-based probability perturbation method (PPM) is applied to characterize uncertain discrete fracture network parameters. The posterior probability distributions of primary and secondary fracture transmissivity, aperture, length, height, and intensity of the secondary fracture are estimated according to the production (flow) histories. A case study based on the Horn River shale gas reservoir is presented. A pilot point scheme is formulated to update the distributions of P32L. The forward modelling entails upscaling each fracture model into a dual-porosity dual-permeability model and performing multiphase flow simulation. This approach produces an ensemble of history-matched discrete fracture and upscaled models by incorporating static and dynamic data. The results demonstrate the utility of the developed method for estimating secondary fracture parameters, which are not inferrable from other static information alone. The inference of both primary and secondary fractures offers a more comprehensive characterization of the fracture networks in tight or shale gas reservoirs. Integrating a pilot point parameterization technique and a sequential simulation approach into the PPM framework presents significant novelties in contributing to a comprehensive and effective method of simulating fluid flow in complex fracture systems, addressing the challenges associated with non-linear relationships between fracture model parameters and the corresponding flow response.