The nearby, massive, runaway star Zeta Ophiuchi has a large bow shock detected in optical and infrared, and, uniquely among runaway O stars, diffuse X-ray emission is detected from the shocked stellar wind. Here we make the first detailed computational investigation of the bow shock of Zeta Ophiuchi, to test whether a simple model of the bow shock can explain the observed nebula, and to compare the detected X-ray emission with simulated emission maps. We re-analysed archival {\it Chandra} observations of the thermal diffuse X-ray emission from the shocked wind region of the bow shock, finding total unabsorbed X-ray flux (0.3-2 keV band) corresponding to a diffuse luminosity of \(L_\mathrm{X}=2.33~(0.79-3.45)\times10^{29}\)ergs\(^{-1}\). 3D MHD simulations were used to model the interaction of the star's wind with a uniform ISM using a range of stellar and ISM parameters motivated by observational constraints. Synthetic infrared, Ha, soft X-ray, emission measure, and radio 6\,GHz emission maps were generated from three simulations, for comparison with relevant observations. Simulations where the space velocity of Zeta Ophiuchi has a significant radial velocity produce infrared emission maps with opening angle of the bow shock in better agreement with observations than for the case where motion is fully in the plane of the sky. The simulation with the highest pressure has the closest match, with flux level within a factor of 2 of the observational lower limit, and emission weighted temperature of \(\log_{10}(T_\mathrm{A}/\mathrm{K})=6.4\), although the morphology of the diffuse emission appears somewhat different. Observed X-ray emission is a filled bubble brightest near the star whereas simulations predict brightening towards the contact discontinuity as density increases.