<p>The study of the flow in a single fracture is the starting point to understand the complex hydraulic behaviour of geological formations and fractured reservoirs, whose comprehension is of interest in many natural phenomena (e.g., magma intrusion) and the optimization of numerous industrial activities in fractured reservoirs (e.g., Enhanced Oil Recovery, drilling engineering, geothermal energy exploitation). Despite the considerable technical prospects of this topic, the associated mathematical complexity and computational burden have so far mostly discouraged investigations of the combined effects of fracture heterogeneity and of the complex rheology of relevant fluids. Indeed, magmas, foams, muds, and suspensions of natural colloids such as clay particles in water are complex fluids and often present in subsurface applications and natural processes. These fluids are characterized by a shear-thinning behavior, which can be well described by the Ellis model, a continuous three-parameter model that behaves as a power-law fluid at high shear rates and as a Newtonian fluid at low shear rates. The Ellis model parameters are: <em>n</em> the power law exponent, <em>&#956;</em><sub>0</sub> the low shear rates viscosity, and <em>&#964;</em><sub>1/2</sub> the shear rate such that <em>&#956;<sub>app</sub></em>(<em>&#964;</em><sub>1/2</sub>)=<em>&#956;</em><sub>0</sub>/2. We use this rheological description in combination with the lubrication theory, which is a depth-averaged formalism permitting us to reduce the full 3-D problem to a 2-D plane formulation. It has been applied to study Newtonian flow in a single fracture for decades and, as far as the aperture gradient remains small (&#8711;<em>d</em>&#171;1), the approximation error introduced by this model is limited. We present here a lubrication-based numerical code aiming at simulating the flow of an Ellis fluid in rough-walled fractures. The code is composed of two modules: a 2D FFT-based fracture aperture field generator and a lubrication-based non-Newtonian flow solver. The former module generates a random aperture field <em>d</em>(<em>x</em>,<em>y</em>) with isotropic spatial correlations, given a mean aperture &#10216;<em>d</em>&#10217;, a coefficient of variation <em>&#963;<sub>d</sub></em>/&#10216;<em>d</em>&#10217;, a Hurst exponent (<em>H</em>) and a correlation length (<em>l<sub>c</sub></em>), reproducing realistic geometries of geological fractures. In the latter module, a 2-D finite volume scheme is adopted to solve the non-linear lubrication equation describing the flow of an Ellis fluid. The equation is discretized on a staggered grid, so that <em>d</em>(<em>x</em>,<em>y</em>) and the pressure field <em>p</em>(<em>x</em>,<em>y</em>)&#160;are defined at different locations. Computational efficiency is achieved by means of the inexact Newton algorithm, with the linearized symmetric system of equations solved via variable-fill-in Incomplete Cholesky Preconditioned Conjugate Gradient method (ICPCG), and a parameter-continuation strategy for the cases with strong nonlinearities. The code proves to be stable and robust when solving flow within strongly heterogeneous fractures (e.g., <em>&#963;<sub>d</sub></em>/&#10216;<em>d</em>&#10217;=1), even on very fine and coarse meshes (e.g., 2<sup>14</sup>&#215;2<sup>14</sup>) and considering a wide range of power-law exponents (e.g., 0.1&#8804;<em>n</em>&#8804;1). The code is validated by comparing the results against analytical solutions (e.g., parallel plates model, sinusoidal profile) and full 3-D CFD simulations, considering different closures.</p>