Snow is an important climate regulator because it greatly increases the surface albedo of middle and high latitudes of the Earth. Earth system models (ESMs) often adopt two-stream approximations with different radiative transfer techniques, the same snow therefore has different solar radiative properties depending whether it is on land or on sea ice. Here we intercompare three two-stream algorithms widely used in snow models, improve their predictions at large zenith angles, and introduce a hybrid model suitable for all cryospheric surfaces in ESMs. The algorithms are those employed by the SNow ICe and Aerosol Radiative (SNICAR) module used in land models, dEdd–AD used in Icepack, the column physics used in the Los Alamos sea ice model CICE and MPAS-Seaice, and a two-stream discrete-ordinate (2SD) model. Compared with a 16-stream benchmark model, the errors in snow visible albedo for a direct-incident beam from all three two-stream models are small (<±0.005) and increase as snow shallows, especially for aged snow. The errors in direct near-infrared (near-IR) albedo are small (<±0.005) for solar zenith angles θ<75∘, and increase as θ increases. For diffuse incidence under cloudy skies, dEdd–AD produces the most accurate snow albedo for both visible and near-IR (<±0.0002) with the lowest underestimate (−0.01) for melting thin snow. SNICAR performs similarly to dEdd–AD for visible albedos, with a slightly larger underestimate (−0.02), while it overestimates the near-IR albedo by an order of magnitude more (up to 0.04). 2SD overestimates both visible and near-IR albedo by up to 0.03. We develop a new parameterization that adjusts the underestimated direct near-IR albedo and overestimated direct near-IR heating persistent across all two-stream models for θ>75∘. These results are incorporated in a hybrid model SNICAR-AD, which can now serve as a unified solar radiative transfer model for snow in ESM land, land ice, and sea ice components.