This is the first of three linked papers which develop an eddy parameterization scheme for mean flows which are wide compared with a deformation radius. The scheme is partly based on the behavior of potential vorticity and thickness fluxes in linear instability, where the former are downgradient (apart from a turning matrix, not present in channel models) and the latter are not precisely downgradient, except on an f-plane. The scheme leads to a diffusivity which varies quite strongly with depth and is smallest at surface and floor. Intrinsic delta-function fluxes also occur at surface and floor, and these are worked out in detail. It is shown that all such parameterization schemes (whether linked to linear instability or not) must satisfy a necessary consistency condition, in the form of a vertical integral. A uniform diffusivity does not satisfy this requirement unless it is defined to vanish at surface and floor. Two methods to compute approximate diffusivities efficiently are given, and their results compare well with exact results from instability theory.