Large eddy simulations were performed for a confined rectangular co-flowing liquid jet at Reynolds number 20,000 based on the average velocity and hydraulic diameter of the channel. An incompressible finite-difference formulation of the filtered Navier-Stokes and mass conservation equations on a partially-staggered grid was used. The effects of grid resolution, numerical schemes, and subgrid models on the LES solutions were studied. Validation was performed by comparing LES statistics with those obtained from particle image velocimetry and laser-induced fluorescence measurements. One-point LES statistics such as the mean, velocity and scalar variances, Reynolds stress, correlation coefficients, skewness, and kurtosis compare well with experimental data. Two point spatial correlations that provide structural information were also computed: very good agreement is observed, showing that structures in the LES field are similar to those in the actual flow.