A growing body of applied mathematics literature in recent years has focussed on the application of fractional calculus to problems of anomalous transport. In these analyses, the anomalous transport (of charge, tracers, fluid, etc.) is presumed attributable to long-range correlations of material properties within an inherently complex, and in some cases self-similar, conducting medium. Rather than considering an exquisitely discretized (and computationally intractable) representation of the medium, the complex and spatially correlated heterogeneity is represented through reformulation of the PDE governing the relevant transport physics such that its coefficients are, instead, smooth but paired with fractional-order space derivatives. Here we apply these concepts to the scalar Helmholtz equation and its use in electromagnetic interrogation of Earth's interior through the magnetotelluric method. We outline a practical algorithm for solving the Helmholtz equation using spectral methods coupled with finite element discretization. Execution of this algorithm for the magnetotelluric problem reveals several interesting features observable in field data: long--range correlation of the predicted electromagnetic fields; a power-law relationship between the squared impedance amplitude and squared wavenumber whose slope is a function of the fractional exponent within the governing Helmholtz equation; and, a non-constant apparent resistivity spectrum whose variability arises solely from the fractional exponent. In geologic settings characterized by self--similarity (e.g. fracture systems; thick and richly-textured sedimentary sequences, etc.) we posit that diagnostics are useful for geologic characterization of features far below the typical resolution limit of electromagnetic methods in geophysics.