  Radio synchrotron polarization maps of the Galaxy can be used to infer the
properties of interstellar turbulence in the diffuse warm ionized medium (WIM).
In this paper, we investigate the spatial gradient of linearly polarized
synchrotron emission ($|\nabla\textbf{P}|/|\textbf{P}|$) as a tracer of
turbulence, the relationship of the gradient to the sonic Mach number of the
WIM, and changes in morphology of the gradient as a function of Galactic
position in the southern sky. We use data from the S-band Polarization All Sky
Survey (S-PASS) to image the spatial gradient of the linearly polarized
synchrotron emission ($|\nabla \textbf{P}|/|\textbf{P}|$) of the entire
southern sky at $2.3$~GHz. The spatial gradient of linear polarization reveals
rapid changes of the density and magnetic fluctuations in the WIM due to
magnetic turbulence as a function of Galactic position; we make comparisons of
these data to ideal MHD numerical simulations. In order to constrain the sonic
Mach number ($M_{s}$), we apply a high order moments analysis to the
observations and to the simulated diffuse, isothermal ISM with ideal
magneto-hydrodynamic turbulence. We find that polarization gradient maps reveal
elongated structures, which we associate with turbulence in the ISM. Our
analysis corroborates the view of a turbulent WIM in a transonic regime
$M_{s}\lesssim 2$. Filamentary structures with typical widths down to the
angular resolution are seen and the observed morphologies match closely with
numerical simulations and in some cases H$\alpha$ contours. The $|\nabla
\textbf{P}|/|\textbf{P}|$ intensity is found to be approximately log-normal
distributed. No systematic variations of the sonic Mach number are observed as
a function of Galactic coordinates, which is consistent with turbulence in the
WIM inferred by the analysis of H$\alpha$ data.
