where $\frac{\partial \mathcal{X}^*}{\partial \hat{d}_i}$ can be expanded as follows:

$$ \frac{\partial \mathcal{X}^*}{\partial \hat{m}_1} := \left[ \begin{array}{l} \frac{\partial \mathcal{X}_{11}}{\partial \hat{m}_1} \cos^2(\hat{\theta}) + \frac{\partial \mathcal{X}_{22}}{\partial \hat{m}_1} \sin^2(\hat{\theta}) \left( \frac{\partial \mathcal{X}_{11}}{\partial \hat{m}_1} - \frac{\partial \mathcal{X}_{22}}{\partial \hat{m}_1} \right) \sin(\hat{\theta}) \cos(\hat{\theta}) \\ \left( \frac{\partial \mathcal{X}_{11}}{\partial \hat{m}_1} - \frac{\partial \mathcal{X}_{22}}{\partial \hat{m}_1} \right) \sin(\hat{\theta}) \cos(\hat{\theta}) \frac{\partial \mathcal{X}_{22}}{\partial \hat{m}_1} \cos^2(\hat{\theta}) + \frac{\partial \mathcal{X}_{11}}{\partial \hat{m}_1} \sin^2(\hat{\theta}) \end{array} \right], \quad (A.8a) $$

$$ \frac{\partial \mathcal{X}^*}{\partial \hat{m}_2} := \left[ \begin{array}{l} \frac{\partial \mathcal{X}_{11}}{\partial \hat{m}_2} \cos^2(\hat{\theta}) + \frac{\partial \mathcal{X}_{22}}{\partial \hat{m}_2} \sin^2(\hat{\theta}) \left( \frac{\partial \mathcal{X}_{11}}{\partial \hat{m}_2} - \frac{\partial \mathcal{X}_{22}}{\partial \hat{m}_2} \right) \sin(\hat{\theta}) \cos(\hat{\theta}) \\ \left( \frac{\partial \mathcal{X}_{11}}{\partial \hat{m}_2} - \frac{\partial \mathcal{X}_{22}}{\partial \hat{m}_2} \right) \sin(\hat{\theta}) \cos(\hat{\theta}) \frac{\partial \mathcal{X}_{22}}{\partial \hat{m}_2} \cos^2(\hat{\theta}) + \frac{\partial \mathcal{X}_{11}}{\partial \hat{m}_2} \sin^2(\hat{\theta}) \end{array} \right], \quad (A.8b) $$

$$ \frac{\partial \mathcal{X}^*}{\partial \hat{\theta}} := \left[ \begin{array}{ll} \sin(\hat{\theta}) \cos(\hat{\theta}) (2\mathcal{X}_{22} - 2\mathcal{X}_{11}) & (\mathcal{X}_{22} - \mathcal{X}_{11}) (\sin^2(\hat{\theta}) - \cos^2(\hat{\theta})) \\ (\mathcal{X}_{22} - \mathcal{X}_{11}) (\sin^2(\hat{\theta}) - \cos^2(\hat{\theta})) & \sin(\hat{\theta}) \cos(\hat{\theta}) (2\mathcal{X}_{11} - 2\mathcal{X}_{22}) \end{array} \right]. \quad (A.8c) $$

Likewise, $\frac{\partial \mathcal{K}^*}{\partial \hat{d}_i}$ can be expanded using the same formulation as described above, and will be elaborated on here.

The weak form of the adjoint equations derived above are given here. As for the objective function, the adjoint equations are as follows:
Find $T_a \in \mathcal{P}$ with $T_a = 0$ on $\Gamma_{in}$ such that:

$$ \begin{aligned} & \int_D - (\epsilon^2 \mathcal{X}^* \nabla p \cdot \nabla \delta T) T_a \, d\Omega + \int_D \mathcal{K}^* \nabla \delta T \cdot \nabla T_a \, d\Omega \\ & = - \int_D \hat{m}_1, \hat{m}_2 \delta T \, d\Omega \quad \forall \delta T \in \mathcal{P}, T_a \in \mathcal{P}, \end{aligned} \quad (A.9a) $$

Find $p_a \in \mathcal{P}$ with $p_a = 0$ on $\Gamma_{out}$ such that:

$$ \int_D \mathcal{X}^* \nabla p_a \cdot \nabla \delta p \, d\Omega = \int_D (\epsilon^2 \mathcal{X}^* \nabla \delta p \cdot \nabla T) T_a \, d\Omega \quad \forall \delta p \in \mathcal{P}, p_a \in \mathcal{P}. \quad (A.9b) $$

As for the pressure drop constraint, $T_a = 0$ in $D$. Thus, the adjoint system reads:
Find $p_a \in \mathcal{P}$ with $p_a = 0$ on $\Gamma_{out}$ such that:

$$ \int_D \mathcal{X}^* \nabla p_a \cdot \nabla \delta p \, d\Omega = -\epsilon^2 \frac{1/\mathrm{DP}_0}{|\Gamma_{\mathrm{in}}|} \int_{\Gamma_{\mathrm{in}}} \delta p \, d\Gamma \quad \forall \delta p \in \mathcal{P}, p_a \in \mathcal{P}. \quad (A.10) $$

## A.2 Discrete adjoint method

We denote the residual of the governing equations as $r_{Darcy}$ and $r_{Energy}$ for the Darcy equation and convection diffusion equation. The Lagrange function $\mathcal{L}$ is constructed as follows:

$$ \mathcal{L} := J + \lambda_p^T r_{Darcy} + \lambda_T^T r_{Energy}, \quad (A.11) $$
where $\lambda_p$ and $\lambda_T$ are the vectors of adjoint variables $p_a$ and $T_a$, respectively. The total derivative of the Lagrange function w.r.t. $\hat{d} = (\hat{m}_1, \hat{m}_2, \hat{\theta})$ is then taken of the Lagrange function, as follows:

$$ \frac{d\mathcal{L}}{d\hat{d}} = \frac{dJ}{d\hat{d}} + \lambda_p^T \frac{d\mathbf{r}_{\text{Darcy}}}{d\hat{d}} + \lambda_T^T \frac{d\mathbf{r}_{\text{Energy}}}{d\hat{d}}, \quad (A.12) $$

where the total derivative of the cost function is expressed as follows:

$$ \frac{dJ}{d\hat{d}} = \frac{\partial J}{\partial \hat{d}} + \frac{\partial J}{\partial \mathbf{s}_p} \frac{\partial \mathbf{s}_p}{\partial \hat{d}} + \frac{\partial J}{\partial \mathbf{s}_T} \frac{\partial \mathbf{s}_T}{\partial \hat{d}}, \quad (A.13) $$

due to the implicit dependence of cost function $J$ (and $G_i$) on the state variables $\mathbf{s}$. Thus, the total derivative of the Lagrange function can be expanded as follows:

$$ \begin{aligned} \frac{d\mathcal{L}}{d\hat{d}} = & \frac{\partial J}{\partial \hat{d}} + \frac{\partial J}{\partial \mathbf{s}_p} \frac{\partial \mathbf{s}_p}{\partial \hat{d}} + \frac{\partial J}{\partial \mathbf{s}_T} \frac{\partial \mathbf{s}_T}{\partial \hat{d}} \\ & + \lambda_p^T \left( \frac{\partial \mathbf{r}_{\text{Darcy}}}{\partial \hat{d}} + \frac{\partial \mathbf{r}_{\text{Darcy}}}{\partial \mathbf{s}_p} \frac{\partial \mathbf{s}_p}{\partial \hat{d}} \right) \\ & + \lambda_T^T \left( \frac{\partial \mathbf{r}_{\text{Energy}}}{\partial \hat{d}} + \frac{\partial \mathbf{r}_{\text{Energy}}}{\partial \mathbf{s}_T} \frac{\partial \mathbf{s}_T}{\partial \hat{d}} + \frac{\partial \mathbf{r}_{\text{Energy}}}{\partial \mathbf{s}_p} \frac{\partial \mathbf{s}_p}{\partial \hat{d}} \right), \end{aligned} \quad (A.14) $$

which can be rewritten as follows:

$$ \begin{aligned} \frac{d\mathcal{L}}{d\hat{d}} = & \frac{\partial J}{\partial \hat{d}} + \lambda_p^T \frac{\partial \mathbf{r}_{\text{Darcy}}}{\partial \hat{d}} + \lambda_T^T \frac{\partial \mathbf{r}_{\text{Energy}}}{\partial \hat{d}} \\ & + \left( \frac{\partial J}{\partial \mathbf{s}_p} + \lambda_p^T \frac{\partial \mathbf{r}_{\text{Darcy}}}{\partial \mathbf{s}_p} + \lambda_T^T \frac{\partial \mathbf{r}_{\text{Energy}}}{\partial \mathbf{s}_p} \right) \frac{\partial \mathbf{s}_p}{\partial \hat{d}} \\ & + \left( \frac{\partial J}{\partial \mathbf{s}_T} + \lambda_T^T \frac{\partial \mathbf{r}_{\text{Energy}}}{\partial \mathbf{s}_T} \right) \frac{\partial \mathbf{s}_T}{\partial \hat{d}}. \end{aligned} \quad (A.15) $$

The adjoint equations are then defined as what are inside the brackets, as follows:

$$ \left( \frac{\partial \mathbf{r}_{\text{Energy}}}{\partial \mathbf{s}_T} \right)^{\mathrm{T}} \boldsymbol{\lambda}_T = - \frac{\partial J}{\partial \mathbf{s}_T}, \quad (A.16a) $$

$$ \left( \frac{\partial \mathbf{r}_{\text{Darcy}}}{\partial \mathbf{s}_p} \right)^{\mathrm{T}} \boldsymbol{\lambda}_p = - \frac{\partial J}{\partial \mathbf{s}_p} - \lambda_T^T \frac{\partial \mathbf{r}_{\text{Energy}}}{\partial \mathbf{s}_p}. \quad (A.16b) $$

Finally, the sensitivity w.r.t. $\hat{\gamma}$ can be derived as follows:

$$ \frac{d\mathcal{L}}{d\hat{d}} = \frac{\partial J}{\partial \hat{d}} + \lambda_p^T \frac{\partial \mathbf{r}_{\text{Darcy}}}{\partial \hat{d}} + \lambda_T^T \frac{\partial \mathbf{r}_{\text{Energy}}}{\partial \hat{d}}. \quad (A.17) $$