Quantitative models are developed for representing the global distribution of the average magnetic field produced by the region 1 and 2 Birkeland current systems. The problem is solved in four following steps: (1) constructing a realistic hit-dependent model of the Birkeland current sheets, based on the formalism of Euler potentials, (2) numerically computing their field at a large number of points within the modeling region, (3) finding a best-fit analytical approximation for that field, and (4) adding a current-free shielding held which confines the Birkeland held within the model magnetopause. At low altitudes the model field-aligned currents reach the ionosphere along eccentric ovals, which fit the observed region I and 2 zones of Iijima and Potemra, and they continue there as horizontal currents. At larger distances the nightside region I currents map to the plasma sheet boundary layer and are then diverted toward the tail flanks, while currents in the dawn-dusk and dayside sectors connect directly to the higher-latitude magnetopause. The region 2 current closes azimuthally near the equator, forming a spread-out partial ring current system. The described approach allows a great flexibility in the geometry of the Birkeland currents, making it feasible to infer their properties from spacecraft data.