The theory of surface operators is described and applied to four surface operators on spheres: the dimensionless surface gradient, ▽1=r▽‐r∂r; the dimensionless surface curl,Λ1=rˆ×▽1; the dimensionless surface Laplacian, ▽1²=▽1· ▽1; and the Funk‐Hecke operators, integral operators with axisymmetric kernels. Three methods are given for solving ▽1²g=fasg=▽1−2f; one method works numerically whenfhas a rapidly convergent expansion in spherical harmonics, the second works whenfis smooth in longitude but not latitude, and the third (a Funk‐Hecke operation) works whenfis rough in all directions. With this apparatus, a complete proof is given of the Helmholtz representation of an arbitrary vector fieldvS(r), the spherical surface of radiusrcentered on the origin: there are unique scalar fieldsf, g, honS(r) such thatv=rˆf+▽1g+Λ1hand 〈g〉r=〈h〉r=0. Here 〈g〉ris the average value ofgonS(r). From the Helmholtz representation on spherical surfaces, the Mie or poloidal‐toroidal representation in spherical shells is deduced. SupposeS(a,c) is the spherical shell whose inner and outer boundaries areS(a) andS(c). SupposeBis solenoidal inS(a,c), i.e., ▽·B=0 and 〈Br〉a=0. Then there are unique scalar fieldsPandQinS(a,c) such thatB=▽ × Λ1P+ Λ1Qand 〈P〉r=〈Q〉r=0 fora⪕r⪕c. The fieldsP= ▽ × Λ1PandQ=Λ1Qare the poloidal and toroidal parts ofB. Applications of this formalism to geomagnetic field modeling are discussed. Gauss's resolution of the geomagnetic fieldBonS(b) into internal and external parts is generalized; if the radial currentJrdoes not vanish onS(b), then to Gauss's expression must be added a toroidal field onS(b) due entirely toJronS(b). A simple proof is given of Runcorn's theorem that to first order in susceptibility no external magnetic field results from magnetization in a horizontally homogeneous spherical shell polarized by sources inside the shell. A Funk‐Hecke‐based method of modeling ionospheric currents is described, which may be more accurate than truncated spherical harmonic expansions and easier to use than Biot‐Savart integrals. Finally, the formalism makes possible the modeling of satellite samples of the geomagnetic field in a spherical shellS(a,c) where electric currents cannot be neglected. Two approximation schemes are described. One is a truncated power series expansion in (c‐a)/H, whereHis the radial length scale of the currents. The other assumes that most ofBinS(a,c) is not due to the currents betweenS(a) andS(c), and that the currents inS(a,c) are field‐aligned. Then the collection of physically possible magnetic fields inS(a,c) is only 50% larger, in a well‐defined sens