The existence of vugs or cavities in naturally fractured reservoirs has long been observed. Even though these vugs can be largely attributed to reserves of oil, natural gas, and groundwater, few investigations of vuggy fractured reservoirs have been conducted. In this paper, a new multi-continuum conceptual model is developed, based on geological data and observations of core examples from carbonate formations in China, to investigate multiphase flow behavior in such vuggy fractured reservoirs. The conceptual model has been implemented into a three-dimensional, three-phase reservoir simulator with a generalized multi-continuum modeling approach. The conceptual model considers vuggy fractured rock as a triple-continuum medium, consisting of (1) highly permeable fractures, (2) low-permeability rock matrix, and (3) various-sized vugs. The matrix system may contain a large number of small or isolated cavities (of centimeters or millimeters in diameter), whereas vugs are larger cavities, with sizes from centimeters to meters in diameter, indirectly connected to fractures through small fractures or microfractures. Similar to the conventional double-porosity model, the fracture continuum is responsible for the occurrence of global flow, while vuggy and matrix continua, providing storage space, are locally connected to each other (and interacting with globally connecting fractures). For practical application of the multi-continuum concept in reservoir simulation, we propose a novel upscaling method for computing equivalent gridblock permeabilities of coarse blocks containing large isolated vugs, in which the local problems consisting of Darcy and Stokes flows are solved. In addition, we describe an efficient boundary condition for accurate computation of upscaled permeabilities. In the numerical implementation, a control-volume, integral finite-difference method is used for spatial discretization, and a first-order finite-difference scheme is adapted for temporal discretization of governing flow equations in each continuum. The resulting discrete nonlinear equations are solved fully implicitly by Newton iteration. The numerical scheme is verified and applied to simulate water-oil flow through the fractured vuggy reservoirs of Tahe Oil Field in China.