Boundary Term in $f(R,T)$ Gravity

Question

In standard $f(R)$ gravity we consider the Lagrangian of the form $L=\frac{1}{16\pi G}f(R)\epsilon$, where $\epsilon$ is the spacetime volume form and similarly, we consider the boundary term to be of the form $l=\frac{1}{8\pi G}f'(R)\epsilon_{\partial \mathcal{M}}$, where $\epsilon_{\partial \mathcal{M}}$ is the spacetime boundary form. Now, in $f(R,T)$ consider the following action
$$S = \int_{\mathcal{M}}L +\int_{\partial \mathcal{M}}l$$
where $L= \frac{1}{16\pi G}f(R,T)\epsilon$ and $l=\frac{1}{8\pi G}f'(R,T)\epsilon_{\partial \mathcal{M}}$, where $f'(R,T)=f_{R}\delta R+f_{T}\frac{\delta\left(g_{\alpha \beta}T^{\alpha \beta}\right)}{\delta g_{\mu \nu}}\delta g_{\mu \nu}$. Upon varying the Lagrangian we would obtain the following form

$$\delta L =\underbrace{\frac{1}{16 \pi G}\left(-R^{\mu\nu}f_{R}+f_{T}\frac{\delta\left(g_{\alpha \beta}T^{\alpha \beta}\right)}{\delta g_{\mu \nu}}\delta g_{\mu \nu} + \frac{1}{2}g^{\mu\nu}f\right)\epsilon\cdot \delta g_{\mu \nu}}_{=E^{\mu\nu}\delta g_{\mu \nu}} + d\Theta,$$

where

$$\theta^{\mu} = \frac{1}{16\pi G}\left(g^{\mu\nu}\nabla^{\nu}g_{\alpha \nu} - g^{\alpha\beta}\nabla^{\mu}\delta g_{\alpha \beta}\right)f_{R},$$

such that $\Theta = \theta\cdot \epsilon$. Now, the variation of the boundary term $l$ is quite messy and considering $\delta f'(R, T) = f_{RR}\delta R +f_{TT}\delta T +f_{RT}\left(\delta R + \delta T \right)$, terms with $\left(\delta g_{\mu\nu}\right)^{2}$ appear which cause problems in trying to fix the pullback of the variation of the metric tensor to the spatial slice after decomposing the boundary term. Computing $\Theta|_{\partial \mathcal{M}}+\delta l$ and imposing Dirichlet's boundary condition, i.e., fixing the pullback of $g_{\mu\nu}$ to $\Gamma$, the stationarity requirement $\left(\Theta +\delta l \right)|_{\Gamma} = dC$, where $C$ is a local $(d-2)$-form on $\Gamma$, is to be fixed. Thus, finally we have the variation of the action to be
$$\delta S = \int_{\mathcal{M}}E^{\mu\nu}\delta g_{\mu\nu} +\int_{\Sigma_{+}\Sigma_{-}}\left(\Theta +\delta l -dC\right), $$

where the following decomposition has been done: $\partial \mathcal{M} = \Gamma\cup\Sigma_{-}\cup\Sigma_{+}$. Firstly, is my boundary term correct, and how am I to fix the squared metric tensor variation term.