It is possible to input a distortion field to 3dNwarpApply, computed from a field map, to transform the EPI datasets. I have tried this, and gotten semi-reasonable results.

However, I'm more pleased with results computed from using reversed blip direction EPI scans, as is done in the Human Connectome Project. There, the 2 sets of EPI datasets are distorted in opposite directions, and then a "meet-in-the-middle" warp can be computed, which will bring them into alignment together, and that works pretty well. But of course that requires changing the way you acquire datasets -- you don't need perfectly balanced acquisitions as the HCP does (half of them blip-left and half of them blip-right), but you do need at least one reference scan of the opposite polarity.