Summary: The geological complexity of fractured reservoirs requires the use of simplified models for flow simulation. This is often addressed in practice by using flow modeling procedures based on the dual-porosity, dual-permeability concept. However, in most existing approaches, there is not a systematic and quantitative link between the underlying geological model [in this case, a discrete fracture model (DFM)] and the parameters appearing in the flow model. In this work, a systematic upscaling procedure is presented to construct a dual-porosity, dual-permeability model from detailed discrete fracture characterizations. The technique, referred to as a multiple subregion (MSR) model, represents an extension of an earlier method that did not account for gravitational effects. The subregions (or subgrid) are constructed for each coarse block using the iso-pressure curves obtained from local pressure solutions of a discrete fracture model over the block. The subregions thus account for the fracture distribution and can represent accurately the matrix-matrix and matrix-fracture transfer. The matrix subregions are connected to matrices in vertically adjacent blocks (as in a dual-permeability model) to capture phase segregation caused by gravity. Two-block problems are solved to provide fracture-fracture flow effects. All connections in the coarse-scale model are characterized in terms of upscaled transmissibilities, and the resulting coarse model can be used with any connectivity-based reservoir simulator. The method is applied to simulate 2D and 3D fracture models, with viscous, gravitational, and capillary pressure effects, and is shown to provide results in close agreement with the underlying DFM. Speedups of approximately a factor of 120 are observed for a complex 3D example.