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 Σ3 109.5° 011 {111} twin boundaries in α-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. cient statistics over a large number of GBs. Consequently, the study of GB segregation usually resorts to a combination of experimental techniques [7, 8] and results depend on the annealing temperature, bulk chemical composition, chemical environment and GB structure. The task becomes even more complicated when studying cosegregation effects, i.e. the simultaneous segregation of at least two chemical species which occurs all the time in real-life materials. There are several mechanisms at work: solute-GB interaction, solute-solute interaction at the GB and in the bulk, site competition, concentration ratio between solutes and kinetic properties of each solute [9]. For instance in Fe-P-C alloys, it was noted that intergranular P segregation decreases with increasing C content but the reason why this happens is not clear [8,10]. If Cr is added to this alloy it seems that P segregation increases only when C concentration is above a certain level [10].The pioneering work of Wu et al.[6] on Σ3[110](111) GB in α-Fe showed that ab initio methods were able to provide some insight on solute properties at GBs, as the authors were able to explain the strengthening by B additions and the embrittlement by P additions in terms of electronic structure at the GB and at free surfaces. These methods are fast-growing, and computation capabilities have increased a lot such that ab initio studies on GBs