Yes—your congruence is an immediate consequence of the classical congruence for Dedekind sums 12 k s(h,k) ≡ h + h^{-1} (mod k), valid for coprime h,k, which itself follows directly from the reciprocity law. You can find this (with your convention for s) in Rademacher–Grosswald, Dedekind Sums, or derive it in two lines from the reciprocity formula s(h,k) + s(k,h) = −1/4 + (1/12)(h/k + k/h + 1/(hk)) by reducing modulo k.

Here is a short derivation of your n=3 statement from that congruence.

Set a' = gcd(b,c), b' = gcd(a,c), c' = gcd(a,b) and d = a' b' c'. Then gcd(bc/d, a/(b'c')) = 1 (and cyclically), so we may apply the basic congruence to each term of S(a,b,c) = a' s(bc/d, a/(b'c')) + b' s(ac/d, b/(a'c')) + c' s(ab/d, c/(a'b')).

For the first term we get 12·(a/(b'c'))·a'·s(bc/d, a/(b'c')) ≡ a'[ (bc/d) + (bc/d)^{-1} ] (mod a'·a/(b'c')). Multiplying by abc / (a/(b'c')) = abc·(b'c')/a and reducing modulo abc gives 12 abc·a' s(bc/d, a/(b'c')) ≡ (bc)^2 + bc·d·(bc/d)^{-1} (mod abc). Summing the three cyclic congruences yields 12 abc S(a,b,c) ≡ (ab)^2 + (bc)^2 + (ca)^2 + d·T (mod abc), where T = bc·(bc/d)^{-1} + ca·(ca/d)^{-1} + ab·(ab/d)^{-1}, and each inverse is taken modulo a/(b'c'), b/(a'c'), c/(a'b') respectively.

Now reduce T modulo a, b, c. For example, from (bc/d)·(bc/d)^{-1} ≡ 1 (mod a/(b'c')) we obtain bc·(bc/d)^{-1} ≡ d (mod a), while the other two terms are multiples of a. Hence T ≡ d (mod a), and similarly T ≡ d (mod b) and T ≡ d (mod c). Since gcd(a,b,c) = 1, this implies T ≡ d (mod lcm(a,b,c)) = abc/d, and therefore d·T ≡ d^2 (mod abc). This proves 12 abc S(a,b,c) ≡ (ab)^2 + (bc)^2 + (ca)^2 + d^2 (mod abc).

Thus your congruence is a straightforward corollary of the standard congruence 12 k s(h,k) ≡ h + h^{-1} (mod k).

Generalization to pairwise coprime a1,…,an Let A = a1⋯an and, for each i, set Hi = A/ai = ∏{j≠i} a_j. When the a_i are pairwise coprime, your symmetric sum is S(a1,…,an) = Σ_i s(Hi, ai). Applying 12 ai s(Hi, ai) ≡ Hi + Hi^{-1} (mod ai) and multiplying by A/ai gives 12 A s(Hi, ai) ≡ (A/ai)(Hi + Hi^{-1}) ≡ Hi^2 + (A/ai)·Hi^{-1} (mod A). Summing over i, 12 A S(a1,…,an) ≡ Σ_i Hi^2 + Σ_i (A/ai)·Hi^{-1} (mod A). The second sum is 1 modulo A by the Chinese remainder theorem: for each i, (A/ai) is 0 modulo aj for j≠i and invertible modulo ai, and choosing Hi^{-1} as its inverse modulo ai makes Σ_i (A/ai)·Hi^{-1} ≡ 1 (mod ai) for every i, hence modulo A. Therefore 12 A S(a1,…,an) ≡ Σ_i (∏{j≠i} a_j)^2 + 1 (mod A), which proves your conjecture for n = 4 and, in fact, for all n ≥ 2 in the pairwise coprime case.

Remarks and references

• The key input is the classical congruence 12 k s(h,k) ≡ h + h^{-1} (mod k). See, for example: • H. Rademacher and E. Grosswald, Dedekind Sums, Theorem derived from the reciprocity law. • T. M. Apostol, Modular Functions and Dirichlet Series in Number Theory, discussion of Dedekind sums and their congruences. • K. Girstmair, various papers on congruences of Dedekind sums (e.g., J. Number Theory 117 (2006), 219–229) survey the arithmetic of 12 k s(h,k) modulo k.
• I do not know a place where your particular symmetric 3-variable (or n-variable) combinations are explicitly recorded, but as shown they follow immediately from the one-variable congruence and CRT.