A famous theorem of Selberg asserts that $\log|\zeta(\tfrac12+it)|$ is approximately a normal distribution with mean $0$ and variance $\tfrac12\log\log T$, when we sample $t\in [T,2T]$ uniformly. This extends in a natural way to a plethora of other $L$-functions, one of them being Dirichlet $L$-functions $L(s,\chi)$ with $\chi$ a primitive Dirichlet character. Viewing $\zeta(\tfrac12+it)$ and $L(\tfrac12+it,\chi)$ as normal variables, we expect indepedence between them, meaning that for fixed $V_1,V_2 \in \mathbb{R}$: $$\textrm{meas}_{t \in [T,2T]} \left\{\frac{\log|\zeta(\tfrac12+it)|}{\sqrt{\tfrac12 \log\log T}}\geq V_1 \text{   and   } \frac{\log|L(\tfrac12+it,\chi)|}{\sqrt{\tfrac12 \log\log T}}\geq V_2\right\} \sim \prod_{j=1}^2 \int_{V_j}^\infty e^{-x^2/2} \frac{\textrm{d}x}{\sqrt{2\pi}}.$$
    When $V_j\asymp \sqrt{\log\log T}$, i.e. we are considering values of order of the variance, the asymptotic above breaks down, but the Gaussian behaviour is still believed to hold to order. For such $V_j$ the behaviour of the joint distribution is decided by the moments $$I_{k,\ell}(T)=\int_T^{2T} |\zeta(\tfrac12+it)|^{2k}|L(\tfrac12+it,\chi)|^{2\ell}\, dt.$$ We establish that $I_{k,\ell}(T)\asymp T(\log T)^{k^2+\ell^2}$ for $0<k,\ell \leq 1$. The lower bound holds for all $k,\ell >0$. This allows us to decide the order of the joint distribution when $V_j =\alpha_j\sqrt{\log\log T}$ for $\alpha_j \in (0,\sqrt{2}]$. Other corollaries include sharp moment bounds for Dedekind zeta functions of quadratic number fields, and Hurwitz zeta functions with rational parameter.