Since aquifer parameters may have statistical dependence structures that are present across a huge range of scales, the concepts of fractional Brownian motion (fBm) have been used in both analytic and numerical stochastic settings. Most previous models have used isotropic scaling characterized by a single scalar Hurst coefficient. Any real‐world anisotropy has been handled by an elliptical stretching random K field. We define a d‐dimensional extension of fBm in which the fractional‐order integration may take on different orders in the d primary (possibly nonorthogonal) scaling directions, and the degree of connectivity and long‐range dependence is freely assigned via a probability measure on the unit sphere. This approach accounts for the different scaling found in the vertical versus horizontal directions in sedimentary aquifers and allows very general degrees of continuity of K in certain directions. It also allows for the representation of fracture networks in a continuum setting: The eigenvectors of the scaling matrix describe the primary fracture scaling directions, and discrete weights of fractional integration represent fracture continuity that may be limited to a small number of directions. In a numerical experiment, the motion of solutes through 2‐D “operator‐fractional” Gaussian fields depends very much on transverse Hurst coefficients. Transverse scaling in the range of fractional Gaussian noise engenders greater plume mixing and a transition to a Fickian regime. Higher orders of integration, in the range of fractional Brownian motion, are associated with thicker layering of aquifer sediments and more preferential, unmixed transport. Therefore direct representation of the unique directional scaling properties of an aquifer is important for realistic simulation of transport.