We propose a new statistical physics model to study equilibrium solute segregation at grain boundaries and the resulting embrittlement effect. This low-temperature expansion model is general and efficient, and its parameters can be obtained from atomistic calculations. It is possible to take into account multiple species, multiple segregation sites with different segregation free energies, account for configurational entropy, grain radius and site competition between solutes. As an example, the model is then applied to the study of phosphorus and hydrogen co-segregation at Sigma 3 109. 5 degrees [0 (1) over bar1]{111} twin boundaries in alpha-Fe, using energetic parameters from density-functional theory calculations. We show that P-H interactions may lead to increased P segregation at grain boundaries and cause additional embrittlement compared to the case where P and H are considered separately.