Naturally fractured reservoirs are generally extremely complex. The aim of characterisation and modelling of such reservoirs is to construct numerical models of rock and fractures, preparing input data for reliable stimulation and fluid flow simulation analyses. This requires the knowledge of different fracture heterogeneities and their correlations at well locations and inter-well regions. This study addresses the issues of how to integrate different information from various field data sources and construct comprehensive discrete fracture networks for naturally fractured reservoir. The methodology combines several mathematical and artificial intelligent techniques, which include statistics, geostatistics, fuzzy neural network, stochastic simulation and simulated annealing global optimisation. The study has contributed to knowledge in characterisation and modelling of naturally fractured reservoirs in several ways. It has developed: .An effective and data-dependant fracture characterisation procedure. It examines all the conventional reservoir data sources and their roles towards characterisation of different fracture properties. The procedure has the advantage of being both comprehensive and flexible. It is able to integrate all multi-scaled and diverse fracture information from the different data sources. .An improved hybrid stochastic generation algorithm for modelling discrete fracture networks. The stochastic simulation is able to utilise both discrete and continuum fracture information. It could simulate not only complicated distributions for fracture properties (e.g. multimodal circular statistics and non-parametric distributions) but also their correlations. In addition, with the incorporation of artificial fuzzy neural simulation, discrete multifractal geometry of fracture size and fracture density distribution map could be evaluated and modelled. Compared to most of the previous fracture modelling approach, this model is more flexible and comprehensive. .An improved conditional global optimisation model for modelling discrete fracture networks. The hybrid model takes full advantages of the advanced fracture characterisation using geostatistical and fuzzy neural analyses. Discrete fractures are treated individually and yet continuum information could be modelled. Compared to the stochastic simulation approach, this model produces more representative fracture networks. Compared to the conventional optimisation programs, this model is more versatile and contains superior objective function.