*how significant are these time-boundary effects which warrant $3d$ rather than $d$
You should just check this. Do the simulations with 1 round, 2 rounds, $d$ rounds, $3d$ rounds, and $10d$ rounds. See if you can see jumps in the extrapolated per-round or per-$d$-rounds error rate.
We never really checked if $d$ was sufficient, we just intuitively assumed $3d$ was probably enough. You can do better than us by just actually checking! Probably $3d$ is more than required; inefficient. Anyways, the underlying intuition we had was that, if the computation is conceived as being made up of $d \times d \times d$ blocks, then the $3d$ case gives you three blocks that can actually appear in real computations whereas the $d$ case gives you one fully enclosed block that would never be used in a real computation.
Note that, experimentally, the boundaries really do stand out compared to the bulk. If you look at figure 2 of "Suppressing quantum errors by scaling a surface code logical qubit
" you see two notable features: the spatial boundaries have lower detection fraction (that's the top-vs-bottom group distinction) and the time boundaries have lower detection fraction (that's the dips at the start and the end).
when using $3d$ rounds of stabilizer measurements how do we get the logical error rate of just that middle code cell?
Define $a \oplus b = a(1-b) + b(1-a)$ to be "Bernoulli xoring". If you have a Bernoulli distribution that returns 1 with probability $a$ and a Bernoulli distribution that returns 1 with probability $b$, then you can form a combined Bernoulli distribution that samples each of these sub distributions, adds up their results, and returns 1 if the sum is odd. This combined distribution returns 1 with probability $a \oplus b$.
Consider the equation $p = a \otimes 4 = a \oplus a \oplus a \oplus a$. It's easy to solve for $p$ given $a$, but how do you solve for $a$ given $p$? You could solve the roots of a degree 4 polynomial, or...
Instead of thinking of a Bernoulli distribution as returning 1 with probability $a$, think of it as returning 0 with probability $a^\prime = 1-2a = \overline{2a}$ and returning a 50/50 random result with probability $2a = 1 - a^\prime = \overline{a^\prime}$. This works as long as $a < 50\%$. The benefit of this switch in perspective is that now, when you xor Bernoulli distributions together, the result is random if any one of them triggered their random case and otherwise it is 0. So you can derive that $p^\prime = (a^\prime)^4$. Meaning $\overline{2p} = \overline{2a}^4$ and so $\overline{2p}^{1/4} = \overline{2a}$ giving $a = (1 - (1 - 2p)^{1/4})/2$.
In general, to convert from a repeated-$k$-times Bernoulli into repeated-$m$-times Bernoulli, while preserving the probability-of-xor-being-1, you do
$$p_m = (1 - (1 - 2p_k)^{k/m})/2$$
Note that, if $k \cdot p_k = m \cdot p_m$ is small, this is well approximated by
$$p_m \approx p_k \cdot k/m$$
The method sinter.shot_error_rate_to_piece_error_rate performs the exact conversion, specialized to k=1. Something great about this function is that, because the only place $m$ appears is in that ratio, the 2-to-3 conversion behaves identically to the 1-to-1.5 conversion. For similar reasons, performing a 1-to-4 and then a 1-to-5 conversion is equivalent to a 1-to-20 conversion. $f_b(f_a(p)) = f_{a \cdot b}(p)$.
Beware that if you're not sampling from a Bernoulli distribution, but from an intersection of Bernoulli distributions (e.g. you have 2 observables and they both had to succeed), then all this math changes a bit.
