The flow development of Herschel-Bulkley fluids in a sudden three-dimensional square expansion is studied numerically. The flow is modeled using the mixed-Galerkin finite element formulation to solve the conservation of mass and momentum equations. The Herschel-Bulkley material behavior is described using a regularized model based on the Papanastasiou model. Solutions are obtained for a downstream-to-upstream expansion ratio of 2:1 and for a wide range of pressure gradient values and rheological parameters. The results show that, during the evolution of the flow, two core regions and dead zones at the corners are formed. The extent of the core regions decreases with the pressure gradient and the Reynolds number, and increases with the power-law index. It is also found that the volume flow rate at steady flow increases with the pressure gradient, power-law index, and Reynolds number.