This thesis is concerned with upscaled models for waterflooded naturally fractured reservoirs (NFRs). Naturally fractured petroleum reservoirs provide over 20% of the world’s oil reserves and production. From the fluid-flow point of view, a fractured reservoir is defined as a reservoir in which a number of naturally occurring fractures have a significant effect on reservoir fluid flow. The reservoir rock between the fractures is called the matrix system. Fractured-reservoir simulations completely differ from conventional-reservoir simulations. The challenge of upscaling is to give an accurate representation of the interaction between fractures and matrix blocks. From the geological point of view, fractured reservoirs can exhibit a number of topologically different configurations. These are reservoirs built from (1) matrix blocks that are bounded by fracture planes in all directions (totally fractured reservoirs, TFRs, or sugar cube), (2) matrix blocks that are bounded only by more or less vertical fracture planes (vertically fractured reservoirs, VFRs), and (3) matrix blocks that form a connected domain interdispersed with fractures (partially fractured reservoirs, PFRs). Only the first configuration, which only exceptionally occurs, is considered in conventional simulators. These simulators use the transfer-function and shape-factor approach. The advantage of this approach is that it can quantify, albeit in a semi-empirical way, in a large variety of cases the fracture-matrix interaction. Moreover, it is fast in terms of computational effort. This thesis adopts a more fundamental approach based on an upscaling methodology called homogenization. It has the advantage of allowing a physically more realistic description of recovery from fractured reservoirs. In addition, it can also be more directly related to the geological model, e.g., the configurations mentioned above. We expect that a comparison of the physically more realistic approach (using homogenization) with the conventional (transfer-function) approach will show when the conventional approach is adequate or when the physically more realistic approach should be used. There are both geological and physical constraints for which homogenization or any other upscaling method is impossible. One advantage of homogenization is that it uses a well-defined procedure that is largely based on physical principles to derive upscaled equations. Different upscaled models are obtained when the dimensionless characteristic numbers are of different orders of magnitude with respect to the scaling factor. We have validated the approach by comparing with a fine-grid simulation that incorporates each individual fracture. Chapter 1 gives a critical review of using homogenization and conditions of its validity both from the geological as from the fluid-flow point of view. Chapter 2 introduces homogenization as an upscaling method to obtain the upscaled model for totally fractured reservoirs (TFRs). Moreover, three dimensionless parameters, which characterize the oil production in fractured media, are defined. The most important of these parameters is a Peclet number defined as the ratio between the average time required to imbibe water into the matrix block and the travel time of water in the fracture system. Different fluid-flow regimes based on the Peclet number are defined and quantified. When the Peclet number is low, the rate of imbibition dominates over the convective transport in the fracture. Therefore, the characteristic time of fluid-flow exchange at the boundary between fracture and matrix is short with respect to the residence time of water in the fracture. When the Peclet number is large, the convective transport is dominant over imbibition and hence the oil recovery mechanism is controlled by the rate of imbibition. Also for a large Peclet number, the residence time of water is much shorter than the characteristic time of imbibition. As a result, less fluid-flow exchange occurs for the same amount of injected water. Thus, if we were in the high-Peclet-number regime, the recovery would be low accordingly and the water cut is large. The results of the upscaled model are compared with ECLIPSE and the effective-permeability model results. Chapter 3 deals with the application of homogenization to derive an upscaled model for vertically fractured reservoirs (VFRs). We formulate a fully implicit three-dimensional upscaled numerical model. Furthermore, we develop a computationally efficient numerical approach. The comparison between the fine-grid results and upscaled-VFR results shows that the upscaled-VFR model can give physically reasonable results. As found previously, also here the Peclet number plays an important role. The gravity number plays a secondary role. For low Peclet numbers, the results are sensitive to gravity, but relatively insensitive to the water-injection rate, lateral matrix-column size, and reservoir topology. At high Peclet numbers, the results are relatively insensitive to gravity, but sensitive to the other conditions mentioned above. We conclude that at low Peclet numbers, it is advantageous to increase the water-injection rate to improve the Net Present Value. However, at high Peclet numbers, increasing the flow rate may lead to uneconomical water cuts. Chapter 4 presents a situation where matrix blocks are connected to neighboring blocks so that part of the global flow occurs in the matrix domain. Such a reservoir is also called a partially fractured reservoir (PFR). We apply homogenization to derive an upscaled model for PFRs that combines dual-porosity and dual-permeability concepts simultaneously. We formulate a well-posed fully implicit 3D upscaled numerical model and investigate oil-recovery mechanisms for different dimensionless characteristic numbers. The specific fracture-matrix interface area is introduced as a new dimensionless number. Upscaled-PFR simulations are compared with fine-grid simulations, the TFR model, and ECLIPSE results. For low Peclet numbers and high gravity numbers, the results are sensitive to gravity and water-injection rates, but relatively insensitive to the specific fracture-matrix interface area, matrix-block size, and reservoir topology, i.e., TFR versus PFR. At high Peclet numbers, the results are relatively insensitive to gravity, but sensitive to the other conditions mentioned above. In particular, when the specific fracture-matrix interface area is large, it enhances the imbibition and consequently leads to a higher oil production. If this specific interface area is small, it leads to a considerable retardation of the imbibition process, which leads to an earlier water breakthrough and lower oil recovery. In general, it is shown that for low Peclet numbers, the upscaled-TFR model and the upscaled-VFR model can be replaced without appreciable loss of accuracy by the effective-permeability (single-porosity) model. However, the effective-permeability model shows large deviations from the TFR and VFR model at high Peclet numbers. Moreover, the results show that at low Peclet numbers and high gravity numbers, ECLIPSE (BWR approach) gives poor predictions and overestimates the oil recovery, but at short injection times, shows good agreement with the solution of the TFR, VFR, and PFR model at intermediate Peclet numbers. The ECLIPSE simulations result in inaccurate predictions of the oil production rate at high Peclet numbers. This can be inferred from the discrepancy with respect to the TFR, VFR, and PFR model for which, based on the comparison between the fine-grid model and the upscaled-VFR and upscaled-PFR model, we assert that they accurately predict the oil recovery. Chapter 5 gives a global understanding of flow in fractured reservoirs. We investigate three main physical aspects of multiphase flow in fractured reservoirs, viz., reservoir wettability, viscosity ratio, and heterogeneity in rock-fluid properties. Various aspects of wettability behavior are examined, i.e., the contact angle, mixed wetting, and non-equilibrium effects in capillary pressure. As to heterogeneity, we study the effects of non-uniform matrix-column size, layering, and fracture aperture. We first discuss the results at low Peclet numbers. For only small gravity numbers, the effect of contact angle, delay time for the non-equilibrium-capillary effect, the layer heterogeneity of the matrix-column size, and the matrix permeability can be ignored without appreciable loss of accuracy. The ultimate oil recovery for mixed-wet VFRs is approximately equal to the Amott index, and the oil production does not depend on the absolute value of the phase viscosity, but on its ratio. However, large gravity numbers enhance under-ridding, aggravated by large contact angles, longer delay times, and higher viscosity ratios. Layering can lead to an improvement or deterioration depending on the fracture aperture and permeability distribution. At low Peclet numbers, the fractured reservoir behaves very similar to a conventional (non-fractured) reservoir and depends largely on the viscosity ratio and the gravity number. At high Peclet numbers, after water breakthrough, the oil recovery appears to be approximately proportional to the cosine of the contact angle and inversely proportional to the sum of the oil and water viscosity. In addition, the mixed-wetting effect is more pronounced; there are significant influences of delay time, matrix permeability, matrix-column size, and the column-size distribution on oil recovery. At low gravity numbers and an effective length-to-thickness ratio larger than ten, the oil recovery is independent of the vertical fracture aperture distribution. For the same amount of injected water, the recovery at low Peclet numbers is larger than the recovery at high Peclet numbers. Chapter 6 investigates the idea to reduce the permeability in a fracture that causes short-circuit flow, by injection of a mixture of sulfuric and hydrochloric acid. During the process, reservoir rock dissolves and gypsum subsequently precipitates. We formulate a mathematical and numerical model that simulates co-injection of sulfuric with hydrochloric acid into carbonate rock for both one- and two-dimensional settings. A one-dimensional simulation of the process shows that upstream dissolution is the dominant process, whereas precipitation of gypsum occurs downstream. The two-dimensional simulation shows that downstream fractures are clogged, whereas continuous paths are created upstream. We find that the injection-flow rate and acid ratio are the two main factors influencing precipitation and dissolution patterns within a fracture. The results show qualitative agreement with experimental results obtained in the Dietz-De Josseling laboratory. We conclude that injection of an acid mixture can be used to reduce the short-circuit flow through fractures and thus to improve waterdrive-recovery performance from oil reservoirs that contain fractures between injection and production well. Overall, it is shown that homogenization is a competitive method for upscaling flow in fractured reservoirs. Homogenization can deal with a large variety of geologically realistic fracture geometries. We have given a physical rationale for most of the assumptions that are used in the mathematical literature for the upscaling procedure by homogenization. We have illustrated that if the conditions for upscaling are met, the ensuing model equations are able to deal with all physical aspects of fractured-media flow.