A new technique was developed to estimate the ionospheric total electron content (TEC) from Global Navigation Satellite System (GNSS) satellite signals. The vertically distributed electron density was parameterized by two thin-shell layers (double-shell approach). The spatiotemporal variation of TEC (strictly speaking, partial electron content) associated with each shell was approximated by the functional fitting of spherical surface harmonics. The major improvements over the conventional single-shell approach were as follows: (1) the precise estimation of TEC was achieved; (2) the estimated TEC was less dependent on the choice of shell height; and (3) the equatorial anomaly was captured more correctly. Furthermore, higher and lower shells exhibited a different pattern of local time vs latitude variation, providing information on the ionosphere--thermosphere dynamics.