Characterisation of optimal solutions to second-order Beckmann problem through bimartingale couplings and leaf decompositions