How do the sums ∑ ᵢ₌₁ⁿ∑ ⱼ₌₀ⁿ1/(i+j)² increase and how are they related to the Euler–Mascheroni constant?

The sequence Sn:=∑i=1n∑j=0n1(i+j)2≈γ+1-ln2+lnn+12n-38n2[12pt]{minimal} {amsmath} {wasysym} {amsfonts} {amssymb} {amsbsy} {mathrsfs} {upgreek} {}{-69pt} {document}Sₙ:=∑ ᵢ₌₁ⁿ∑ ⱼ₌₀ⁿ1/(i+j)²≈ γ +1 -ln 2 + ln n+1/2n-3/8n²{document}is estimated sharply using the Euler–Mascheroni constant.