Previous studies have confirmed the adverse impact of fading correlation on the mutual information (MI) of two-dimensional (2D) multiple-input multiple-output (MIMO) systems. More recently, the trend is to enhance the system performance by exploiting the channel's degrees of freedom in the elevation, which necessitates the derivation and characterization of three-dimensional (3D) channels in the presence of spatial correlation. In this paper, an exact closed-form expression for the Spatial Correlation Function (SCF) is derived for 3D MIMO channels. This novel SCF is developed for a uniform linear array of antennas with nonisotropic antenna patterns. The proposed method resorts to the spherical harmonic expansion (SHE) of plane waves and the trigonometric expansion of Legendre and associated Legendre polynomials. The resulting expression depends on the underlying arbitrary angular distributions and antenna patterns through the Fourier Series (FS) coefficients of power azimuth and elevation spectrums. The novelty of the proposed method lies in the SCF being valid for any 3D propagation environment. The developed SCF determines the covariance matrices at the transmitter and the receiver that form the Kronecker channel model. In order to quantify the effects of correlation on the system performance, the information-theoretic deterministic equivalents of the MI for the Kronecker model are utilized in both mono-user and multi-user cases. Numerical results validate the proposed analytical expressions and elucidate the dependence of the system performance on azimuth and elevation angular spreads and antenna patterns. Some useful insights into the behaviour of MI as a function of downtilt angles are provided. The derived model will help evaluate the performance of correlated 3D MIMO channels in the future.