Dual-continuum models are useful for describing flow in porous systems with significant local pressure disequilibrium between slow moving water, contained in the porous matrix, and fast moving water in preferential pathways. The formation and intensity of preferential flow depends on the contrast between the hydraulic properties of the two flow domains as well as on the properties of their interface. In this study, we focused on both physical coupling of the flow domains through the mass transfer term and numerical coupling of the respective governing equations. The set of governing equations was alternatively solved using a sequentially coupled (SC) approach and a fully coupled (FC) approach. The SC approach was shown to be computationally more efficient for strongly developed preferential flow in systems with high interfacial resistance; however, it becomes numerically unstable for weak preferential flow associated with low interfacial resistance. The FC approach is a computationally more expensive yet numerically more robust alternative, capable of simulating a complete class of intermediate flow regimes ranging from strongly preferential flow in a dual-continuum system to nonpreferential flow in a single-continuum system. To illustrate the performance of the numerical coupling approaches in conjunction with the effect of different interfacial resistances, we present a simple example problem involving one-dimensional near-saturated flow in a vertical soil column.