In interest of accurately modelling wellbore production in convective geothermal systems, we draw attention to two approaches to numerical simulation of crustal fluid flow. The historical ''strong'' approach simulates the fluid mass conservation condition on a localised finite difference mesh approximation to poroelastic continua. The more flexible ''weak'' approach globally enforces fluid mass conservation by in effect integrating over all numerical mesh cells. Strong formulation crustal flow simulation emerged from numerical simulations of mathematically tractable conductive heat flow in media typified by normallydistributed/spatially-uncorrelated thermal property distributions. Strong formulation numerical stability problems quickly intrude, however, when confronted with the lognormally-distributed/spatially-correlated crustal permeability distributions that pervade the ambient crust. Weak flow conservation solver formulations accurately handle abrupt wide-ranging crustal permeability spatial fluctuations arising from these crustal permeability