Quantifying modeling uncertainties when combining multiple gravitational-wave detections from binary neutron star sources