A numerical, poro-viscoelastic finite element flow model has been developed to simulate the redistribution of melt during the folding of a horizontal layer. Important mechanical properties governing the overall rheology of the buckling layer, including porosity, matrix viscosity, melt pressure and changes in deviatoric stress are tracked simultaneously with deformation. Melt migrates down gradients in melt pressure to the outer hinge zone of the folds, away from sites of compression to regions of local tension. For a matrix with dynamic power law rheology, deformation results in the formation of low viscosity, strain weakening zones that coincide with zones of enhanced deviatoric stress that are potential sites for plastic deformation. During fast deformation of layers with large width to height ratios, the more rapidly bending material splits into domains of higher viscosity separated by lower viscosity zones. Porosity (melt fraction) distribution patterns do not appear to be controlled by the matrix rheology. The numerical results and predictions of melt distribution and geometry during deformation correspond well with data from some physical (analogue) laboratory experiments and field observations. The quantitative nature of our results should help refine more qualitative models for the rates and mechanisms of melt extraction and transport in folded crustal migmatite terranes.