Analogous to the TV metric, the requested upper bound holds for the discrepancy metric with $k=q$. The result given below can also be easily extended to general product probability measures $\mu=\otimes_{i=1}^q \mu_{i}$ and $\nu=\otimes_{i=1}^q \nu_i$ to obtain $$
D(\mu,\nu)\le \sum_{i=1}^q D(\mu_i, \nu_i)\;.
$$

Suppose that $\mu_1$ and $\mu_2$ are probability measures. Let $B = \{ (x_1, \dots, x_q) : \sum_i (x_i - c_i)^2 \le R^2 \}$ be a closed ball in $\mathbb{R}^q$ with radius $R>0$ centered at $(c_1, \dots, c_q) \in \mathbb{R}^q$. Let $\mathbf{1}_A$ denote the indicator function of the set $A$.

By telescoping, and invoking Tonelli’s theorem to write $(\mu_1^{\otimes q}-\mu_2^{\otimes q})(B)$ as an iterated integral, note that
\begin{align*} & |(\mu_1^{\otimes q}-\mu_2^{\otimes q})(B)|=
\left| \int_{\mathbb{R}^q} \mathbf{1}_B(x_1, \dots, x_q) \biggl( d\mu_1(x_1) \cdots d\mu_1(x_q) - d\mu_2(x_1) \cdots d\mu_2(x_q) \biggr) \right| \\
&= \left| \int_{\mathbb{R}^q} \mathbf{1}_B(x_1, \dots, x_q) \sum_{ i=1}^q d\mu_1(x_1) \cdots d\mu_1(x_{i-1}) \biggl(d\mu_1(x_{ i}) - d\mu_2(x_{ i})\biggr) d\mu_2(x_{i+1}) \dots d\mu_2(x_q) \right| \\
&= \left| \sum_{ i=1}^q \int_{\mathbb{R}^{q-1}} \left(\int_{\mathbb{R}} \mathbf{1}_B(x_1, \dots, x_q)\biggl(d\mu_1(x_{ i}) - d\mu_2(x_{ i})\biggr) \right) d\mu_1(x_1) \cdots d\mu_1(x_{i-1}) d\mu_2(x_{i+1}) \dots d\mu_2(x_q) \right| \\
&\le \sum_{ i=1}^q \int_{\mathbb{R}^{q-1}} \left| \int_{\mathbb{R}} \mathbf{1}_B(x_1, \dots, x_q)\biggl(d\mu_1(x_{ i}) - d\mu_2(x_{ i})\biggr) \right| d\mu_1(x_1) \cdots d\mu_1(x_{i-1}) d\mu_2(x_{i+1}) \dots d\mu_2(x_q) \\
&\le \sum_{i=1}^q \sup_{\gamma_i, \rho_i} \left| \int_{\mathbb{R}}\mathbf{1}_{(x_i-\gamma_i)^2 \le \rho_i^2} \biggl(d\mu_1(x_{ i}) - d\mu_2(x_{ i}) \biggr) \right| \\
& \qquad \int_{\mathbb{R}^{q-1}} d\mu_1(x_1) \cdots d\mu_1(x_{i-1}) d\mu_2(x_{i+1}) \cdots d\mu_2(x_q) \\
&\le q D(\mu_1, \mu_2) \;.
\end{align*}
Since $B$ is arbitrary, the conjectured upper bound holds with $k=q$.

**Remark.** One can do much better than this linear upper bound in some special cases. For instance, suppose that $\mu_1 = \mathcal{N}(0,1)$ and $\mu_2 = \mathcal{N}(0, \sigma^2)$. Then, by moving to hyperspherical coordinates, it's not too hard to show that $$
D(\mu_1^{\otimes q},\mu_2^{\otimes q}) = \left| \frac{\Gamma(\frac{q}{2}, \frac{q \log{\sigma}}{\sigma^2-1}) - \Gamma(\frac{q}{2}, \frac{q \sigma^2 \log(\sigma)}{\sigma^2-1})}{ \Gamma(\frac{q}{2}) } \right| \; $$ which converges to one with $q$ and $k(q):= D(\mu_1^{\otimes q},\mu_2^{\otimes q}) /D(\mu_1,\mu_2) $ grows sublinearly. This Gaussian case is rather exceptional because the corresponding product measure can be directly written in terms of the Euclidean distance.

2more comments