Abstract Natural fracture systems comprise numerous small features and relatively few large ones. At field scale, it is impractical to treat all fractures explicitly. We represent the largest fractures via Embedded Discrete Fracture Modeling (EDFM) and account for smaller ones using a dual-porosity, dual-permeability (DPDK) idealized representation of the fracture network. The hierarchical EDFM+DPDK approach uses consistent discretization schemes and efficiently simulates realistic field cases. Further speed-up can be obtained using aggregation-based upscaling. Capabilities to visualize and post-process simulation results facilitate understanding for effective management of fractured reservoirs. The proposed approach embeds large discrete fractures as EDFM within a DPDK grid (which contains both matrix and idealized fracture continua for smaller fractures), and captures all connections among the triple media. In contrast with existing EDFM formulations, we account for discrete fracture spacing within each matrix cell via a new matrix-fracture transfer term and employ consistent assumptions for classical EDFM and DPDK calculations. In addition, the workflow enables coarse EDFM representations using flow-based cell-aggregation upscaling for computational efficiency, as well as finite-volume tracer-based flux post-processing to analyze production allocation and sweep. Using a synthetic case, we show that the proposed EDFM+DPDK approach provides a close match of simulation results from a reference model that represents all fractures explicitly, while providing runtime speedup. It is also more accurate than previous standard EDFM and DPDK models. We demonstrate that the matrix-fracture transfer function agrees with flow-based upscaling of high-resolution fracture models. Next, the automated workflow is applied to a waterflooding study for a giant carbonate reservoir, with an ensemble of stochastic fracture realizations. The overall workflow provides the computational efficiency needed for performance forecasts in practical field studies, and the 3D visualization allows for the derivation of insights into recovery mechanisms. Finally, we apply a flux post-processing scheme on simulation results to understand expected waterflood performance.