I agree that the $L_2$ norm is off the table, because it creates a nonconvex quadratic constraint. Unfortunately, given the base model (and in particular the choice of binary variables), I'm not sure it would be practical to use the $L_1$ norm either.
The problem with $L_1$ lies in the absolute values, which need to be linearized. There are essentially three ways to handle $\vert u - v\vert$ in a MIP model. The first is easy: if you know that $u > v$, the absolute value is $u-v$, and otherwise it is $v-u$. The problem here is that you don't reliably know the signs of the coordinate differences between modules. For instance, if $x_{ij} = 0 = y_{ij}$, you know that $j$ is to the right of $i$, but you don't know if it is above/alongside/below $i$. The second linearization technique is to use an auxiliary variable $w$ with the constraints $w\ge u - v$ and $w\ge v - u$, which implies $w\ge \vert u-v\vert$. This would be fine if your norm constraint were $\parallel \dots \parallel \le D$, but because it is a greater-than constraint, the model would be tempted to "lie" (e.g., $u=1$, $v=1$, $w=D > 1$). The third approach involves throwing in additional binary variables, one for each absolute value, along with additional defining constraints, which I think would be a huge headache to formulate and would result in a very large model.
A more tractable approach might be to use the $L_\infty$ norm (known in some circles as the "sup" norm). For each pair $i$, $j$ of modules, you need to consider all four possible combinations of $x_{ij}$ and $y_{ij}$. I'll demonstrate one case (the one depicted in your diagram) and leave the others as an exercise for the reader. Suppose we have $x_{ij}=0=y_{ij}$, meaning $i$ is left of $j$, as in the diagram. There are four combinations of a red block from $i$ and a red block from $j$. A sufficient condition that all four are separated by at least $D$ in the $L_\infty$ norm is that the left edge of $j$ is at least $D$ units to the right of the right edge of $i$. Note that this condition is sufficient but not necessary, since $j$ could be closer to $i$ horizontally but, if we slid $j$ up or down enough, further than $D$ from $i$ in the sup norm. To enforce the sufficient condition, we add the constraint $$x_j - (x_i + w_i) \ge D(1-x_{ij}-y_{ij}).$$
Because the approach is "conservative", this might be too wasteful of space, but it is easy to add to the existing model.
I should note that I am assuming all modules have the same layout (blocks 2 and 4 red, 1 and 3 not). If, for instance, module $i$ was as depicted but module $j$ had only blocks 1 and/or 4 red, the left side of the constraint above would become $$\left( x_j + \frac{w_j}{2}\right) - \left( x_i + w_i \right).$$
Update: Sticking to the $L_\infty$ norm, I think there is a way to reduce wasted space (but not eliminate it entirely), at a cost in computational complexity. Whether the trade-off is worthwhile is, to me, an empirical question.
The approach is as follows. First, in addition to continuous variables for the coordinates of the lower left corner of each module, add continuous variables for the lower left corner of each red block, along with linear equations defining the block corners in terms of the module corners. This much is "cheap".
Next, redefine the binary variables to apply to individual red blocks, as well as modules. Since you are going from two binary variables per pair of modules to two per pair of blocks, and there are four times as many block pairs as module pairs, this will quintuple the number of binary variables, as well as the number of distance-related constraints. (The no-overlap constraints can still be enforced at the module level, so no change there.) If $D$ is large enough relative to the module dimensions, you may not need the original binary variables and the no-overlap constraints (because keeping red blocks more than $D$ apart would preclude module overlap), but my guess is you will need to keep the original binaries and their constraints.
