The dominance variance on a single locus
It looks like I could solve the question. My understanding is as follows.

Let $(X_1,Y_1)$, $(X_2,Y_2)$ (corresponds to $i,j$ and $k,l$ in the question) be the genotypes of two individuals 
 at some location. We assume that neither individual
 is in-bred. In this case define,
 \begin{align*}
  \Delta_7&amp;= \text{Pr}(\text{Both alleles $\{X_1,Y_1\}$ and $\{X_2, Y_2\}$ are IBD}),\\
  \Delta_8&amp;= \text{Pr}(\text{Exactly one allele pair of the four pairs is  IBD}),\\
  \Delta_9&amp;= \text{Pr}(\text{ No IBD between the two individuals}).
 \end{align*}
 $\Delta_7$, $\Delta_8$ and $\Delta_9$ are calculated from the pedigrees. 
 For example, for non-twin siblings, $\Delta_7 = \frac{1}{4}$, $\Delta_8=\frac{1}{2}$,
 $\Delta_9 =\frac{1}{4}$.

We note that the IBD coefficient (also called kinship coefficient) $\Theta$ is 
 \begin{align*}
  \Theta = \frac{1}{2} \Delta_7 + \frac{1}{4} \Delta_8.
 \end{align*}
 Also the fraternity coefficient, $\Delta =\Delta_7$.
Thus we have that the covariance between the genotypic 
 value of two individuals is given by
 \begin{align*}
  \mathbb{E}[G(X_1,Y_1)G(X_2,Y_2)] =&amp; \Delta_7 \mathbb{E}[G(X_1,Y_1)G(X_2,Y_2)| \text{both IBD}] + \\
  &amp;\Delta_8 \mathbb{E}[G(X_1,Y_1)G(X_2,Y_2)| \text{one IBD}] + \Delta_9  \mathbb{E}[G(X_1,Y_1)G(X_2,Y_2)| \text{no IBD}].
 \end{align*}

Further note that,
 \begin{align*}
  \mathbb{E}[G(X_1,Y_1)G(X_2,Y_2)| \text{both IBD}] &amp;= \mathbb{E}[G^2(X_1,Y_1)],\\
  &amp;= \sigma_A^2 +\sigma_D^2.
 \end{align*}

\begin{align*}
  \mathbb{E}[G(X_1,Y_1)G(X_2,Y_2)| \text{no IBD}] &amp;= \mathbb{E}[G(X_1,Y_1)] \mathbb{E}[G(X_2,Y_2)],\\
  &amp;=0.
 \end{align*}

\begin{align*}
  \mathbb{E}[G(X_1,Y_1)G(X_2,Y_2)| \text{one IBD}] &amp;= \mathbb{E}[(\alpha(X)+\alpha(Y_1)+\delta(X,Y_1))(\alpha(X)+\alpha(Y_2)+\delta(X,Y_2))],\\
  &amp;=\mathbb{E}[\alpha^2(X)] + 2 \mathbb{E}(\alpha(X) \delta(X,Y_1)) + \mathbb{E} [\delta(X,Y_1)\delta(X,Y_2)],\\
  &amp;= \frac{1}{2} \sigma_A^2 + 2 \mathbb{E} [ \alpha(X) \mathbb{E}[\delta(X,Y_1) | X] ] +  \mathbb{E} [\mathbb{E}[\delta(X,Y_1)|X]\mathbb{E}[\delta(X,Y_2)|X]],\\
  &amp;= \frac{1}{2} \sigma_A^2.
 \end{align*}

Thus we have that,
 \begin{align*}
  \mathbb{E}[G(X_1,Y_1)G(X_2,Y_2)] =&amp; \left(\frac{1}{2}\Delta_7 + \frac{1}{4}\Delta_8\right) 2 \sigma_A^2 + \Delta_7 \sigma_D^2,\\
  &amp;= 2\Theta \sigma_A^2 + \Delta \sigma_D^2,
 \end{align*}
 which is as claimed.