The goal of this study is to determine whether shear wave splitting observed in subduction zone back arc regions, the Tonga subduction zone in particular, can be quantitatively modeled with flow in the back arc mantle driven by the motions of the subducting slab and the upper back arc plate. We calculated two-dimensional mantle flow models using known Tonga plate motions as boundary conditions and assuming a range of uniform and variable viscosity structures. Shear wave splitting was predicted for the anisotropy due to lattice preferred orientation (LPO) of olivine and orthopyroxene in the flow model finite strain fields. The predicted shear wave splitting provides a good match to the fast directions (parallel to the azimuth of subducting plate motion) and splitting times (0.5--1.5 s) observed in Tonga, both for models where LPO anisotropy develops everywhere above 410 km and for models where LPO anisotropy is confined to regions of relatively high stress. If LPO anisotropy does develop over the entire upper 410 km of the mantle, the strength of anisotropy induced by a given amount of shear strain must be relatively weak (~4% for shear strains of 1.5, with a maximum value of ~6% for very large strains). The splitting observations are comparably fit by a wide range of different viscosity models. Anisotropy due to melt-filled cracks aligned by stresses in the back arc flow models predicts fast directions roughly normal to observed values and thus cannot alone explain the observed splitting.