Naturally-fractured carbonate karst reservoirs are characterized by various-sized solution caves that are connected via fracture networks at multiple scales. These complex geologic features can not be fully resolved in reservoir simulations due to the underlying uncertainty in geologic models and the large computational resource requirement. They also bring in multiple flow physics which adds to the modeling difficulties. It is thus necessary to develop a method to accurately represent the effect of caves, fractures and their interconnectivities in coarse-scale simulation models. In this paper, we present a procedure based on our previously proposed Stokes-Brinkman model (SPE 125593) and the discrete fracture network method for accurate and efficient upscaling of naturally fractured carbonate karst reservoirs.