The Sine-Gordon model is obtained by tilting the law of a log-correlated Gaussian field $X$ defined on a subset of $\mathbb{R}^d$ by the exponential of its cosine, namely $\exp(α\smallint \cos (βX))$. It is an important model in quantum field theory or in statistic physics like in the study of log-gases. In spite of its relatively simple definition, the model has a very rich phenomenology. While the integral $\smallint \cos (βX)$ can properly be defined when $β^2<d$ using the standard Wick normalisation of $\cos (βX)$, a more involved renormalization procedure is needed when $β^2\in [d,2d)$. In particular it exhibits a countable sequence of phase transition accumulating to the left of $β=\sqrt{2d}$, each transitions corresponding to the addition of an extra term in the renormalization scheme. The final threshold $β=\sqrt{2}$ corresponds to the Kosterlitz-Thouless (KT) phase transition of the $\log$-gas. In this paper, we present a novel probabilistic approach to renormalization of the two-dimensional boundary (or 1-dimensional) Sine-Gordon model up to the KT threshold $β=\sqrt{2}$. The purpose of this approach is to propose a simple and flexible method to treat this problem which, unlike the existing renormalization group techniques, does not rely on translation invariance for the covariance kernel of $X$ or the reference measure along which $\cos (βX)$ is integrated. To this purpose we establish by induction a general formula for the cumulants of a random variable. We apply this formula to study the cumulants of (approximations of) $\smallint \cos (βX)$. To control all terms produced by the induction proceedure, we prove a refinement of classical electrostatic inequalities, which allows to bound the energy of configurations in terms of the Wasserstein distance between $+$ and $-$ charges.