A logarithmic type modulus of continuity is established for weak solutions to a two-phase Stefan problem, up to the parabolic boundary of a cylindrical space-time domain. For the Dirichlet problem, we merely assume that the spatial domain satisfies a measure density property, and the boundary datum has a logarithmic type modulus of continuity. For the Neumann problem, we assume that the lateral boundary is smooth, and the boundary datum is bounded. The proofs are measure theoretical in nature, exploiting De Giorgi's iteration and refining DiBenedetto's approach. Based on the sharp quantitative estimates, construction of continuous weak (physical) solutions is also indicated. The logarithmic type modulus of continuity has been conjectured to be optimal as a structural property for weak solutions to such partial differential equations.