Can these S-parameters be derived?
I've solved the problem, and it was slightly subtle. It's an interesting problem, so I'll write a detailed solution here for anyone that needs it. This is a lossless transmission-line T-network power divider. 

To answer this question, it's important to note that because the generator and load impedances are, in general, different. If we define the s-matrix from the traveling waves, V+(z) and V-(z), but don't take the different reference impedances into account in these definitions, then our results won't agree with those of generalized scattering parameters. 

In order to obtain results which are consistent, it is necessary to replace V+ and V- variables with values normalized to the reference impedance of that port, often denoted a and b (in the style of Collin). 
$$ a_i=\frac{V_i^+}{\sqrt{Z_{0i}}} $$
$$ b_i=\frac{V_i^-}{\sqrt{Z_{0i}}} $$
The s-parameters are now more properly given by:
$$ S_{ij} = \left. \frac{a_i}{b_j}\right\vert_{V_k=0, k\neq j} = \left.\frac{V_i^+}{V_j^-}\frac{\sqrt{Z_{0j}}}{\sqrt{Z_{0i}}}\right\vert_{V_k=0, k\neq j} $$

Using previously obtained results, it was found that:
$$ V_2^- = -jV_g\frac{Z_{in,p}}{Z_{in,p}+Z_g}\frac{1+\Gamma_{L1}}{1-\Gamma_{L1}} $$
$$ V_1^+ = \frac{1}{2} V_g $$
Which becomes
$$ \frac{V_2^-}{V_1^+} = -j2\frac{Z_{in,p}}{Z_{in,p}+Z_g}\frac{1+\Gamma_{L1}}{1-\Gamma_{L1}} = 2V_1(0) $$
Now normalizing for the different port impedances, we get the value of S21,
$$ S_{21} = 2V_1(0)\sqrt{\frac{Z_g}{Z_{L1}}} $$
Exactly what was obtained empirically in the question.