Different flow situations can be encountered at the boundary between a non-consolidated porous medium and a Newtonian fluid. For a given medium, three flow patterns exist as a function of the flow rate Q, when the fluid flows from the porous medium through the interface: Q< Q 1: filtration through the medium characterised by its permeability. Q 1< Q< Q 2: beyond a local fluidisation yield Q 1, constitutive particles at the interface are torn away and an unstable front develops. Q> Q 2: a second fluidisation yield, global this time, brings about a new homogeneous front corresponding to the movement of the porous medium as a whole. Experiments are made in which we study the flow of a glycerol solution through sand at varying flow rates, in a nearly 2d geometry. As the Q 1 threshold is passed, an erosion instability is quickly established, giving rise to branched waterways penetrating the sand bulk. We propose a theoretical interpretation of the phenomenon, based on a classical perturbation analysis of the water/sand interface. A dispersion equation is found, as well as an equation predicting the number of branches as a function of branch width λ. We find that the waterway patterns obey fractal scaling laws with good precision, as expected in our model, giving a fractal dimension of 1.7 for well-established instabilities.