Abstract Cold heavy oil production with sand (CHOPS) is a non-thermal primary process that is widely adopted in many weakly consolidated heavy oil deposits around the world. However, only 5 to 15% of the initial oil in place is typically recovered. Several solvent-assisted schemes are proposed as follow-up strategies to increase the recovery factor in post-CHOPS operations. The development of complex, heterogeneous, high-permeability channels or wormholes during CHOPS renders the analysis and scalability of these processes challenging. One of the key issues is how to properly estimate the dynamic growth of wormholes during CHOPS. Existing growth models generally offer a simplified representation of the wormhole network, which, in many cases, is denoted as an extended wellbore. Despite it is commonly acknowledged that wormhole growth due to sand failure is likely to follow fractal statistics, there are no established workflows to incorporate geomechanical constraints into the construction of these fractal wormhole patterns. A novel dynamic wormhole growth model is developed to generate a set of realistic fractal wormhole networks during the CHOPS operations. It offers an improvement to the Diffusion Limited Aggregation (DLA) algorithm with a sand-arch-stability criterion. The outcome is a fractal pattern that mimics a realistic wormhole growth path, with sand failure and fluidization being controlled by geomechanical constraints. The fractal pattern is updated dynamically by coupling compositional flow simulation on a locally-refined grid and a stability criterion for the sand arch: the wormhole would continue expanding following the fractal pattern, provided that the pressure gradient at the tip exceeds the limit corresponding to a sand-arch-stability criterion. Important transport mechanisms including foamy oil (non-equilibrium dissolution of gas) and sand failure are integrated. Public field data for several CHOPS fields in Canada is used to examine the results of the dynamic wormhole growth model and flow simulations. For example, sand production history is used to estimate a practical range for the critical pressure gradient representative of the sand-arch-stability criterion. The oil and sand production histories show good agreement with the modeling results. In many CHOPS or post-CHOPS modeling studies, constant wormhole intensity is commonly assigned uniformly throughout the entire domain; as a result, the ensuing models are unlikely to capture the complex heterogeneous distribution of wormholes encountered in realistic reservoir settings. This work, however, proposes a novel model to integrate a set of statistical fractal patterns with realistic geomechanical constraints. The entire workflow has been readily integrated with commercial reservoir simulators, enabling it to be incorporated in practical field-scale operations design.