$\newcommand{\Dev}{\operatorname{Dev}}$ $\newcommand{\dev}{\operatorname{dev}}$ $\newcommand{\grad}{\operatorname{grad}}$ $\newcommand{\Grad}{\operatorname{Grad}}$ $\newcommand{\divg}{\operatorname{div}}$ $\newcommand{\Ei}{\operatorname{Ei}}$ $\newcommand{\tr}{\operatorname{tr}}$ $\newcommand{\Divg}{\operatorname{Div}}$ $\newcommand{\cay}{\operatorname{cay}}$ $\newcommand{\rot}{\operatorname{rot}}$

The contribution to $\delta W$ due to chemical reactions is given by $\delta G$ , where $$\begin{equation} \delta G=\int_{b}\delta\tilde{p}\left[\hat{\varphi}^{w}+\left(1-\varphi^{s}\right)\hat{\zeta}\overline{\mathcal{V}}\right]\,dv+\sum\limits _{\iota}\nu^{\iota}\int_{b}\delta\tilde{c}^{\iota}\left(1-\varphi^{s}\right)\hat{\zeta}\,dv\,.\label{eq354} \end{equation}$$ The linearization of $\delta G$ along a solid displacement increment $\Delta\mathbf{u}$ is $$\begin{equation} \begin{aligned}D\delta G\left[\Delta\mathbf{u}\right] & =\int_{b}\delta\tilde{p}\,\left(\hat{\varphi}^{w}\,\divg\Delta\mathbf{u}+\hat{\boldsymbol{\varphi}}_{\varepsilon}^{w}:\Delta\boldsymbol{\varepsilon}\right)\,dv\\ & +\overline{\mathcal{V}}\int_{b}\delta\tilde{p}\left[\hat{\zeta}\,\divg\Delta\mathbf{u}+\left(J-\varphi_{r}^{s}\right)\hat{\boldsymbol{\zeta}}_{\varepsilon}:\Delta\boldsymbol{\varepsilon}\right]\,dv\\ & +\sum_{\iota}\nu^{\iota}\int_{b}\delta\tilde{c}^{\iota}\left[\hat{\zeta}\,\divg\Delta\mathbf{u}+\left(J-\varphi_{r}^{s}\right)\hat{\boldsymbol{\zeta}}_{\varepsilon}:\Delta\boldsymbol{\varepsilon}\right]\,dv \end{aligned} \,,\label{eq354b} \end{equation}$$ where $$\begin{equation} \hat{\boldsymbol{\varphi}}_{\varepsilon}^{w}=\mathbf{F}\cdot\frac{\partial\hat{\varphi}^{w}}{\partial\mathbf{E}}\cdot\mathbf{F}^{T}\,,\quad\hat{\boldsymbol{\zeta}}_{\varepsilon}=J^{-1}\mathbf{F}\cdot\frac{\partial\hat{\zeta}}{\partial\mathbf{E}}\cdot\mathbf{F}^{T}\,.\label{eq354c} \end{equation}$$ Currently, $\hat{\zeta}$ is assumed to be independent of $\tilde{p}$ in FEBio; it follows that the linearization along the effective fluid pressure increment $\Delta\tilde{p}$ is $$\begin{equation} D\delta G\left[\Delta\tilde{p}\right]=\int_{b}\delta\tilde{p}\,\frac{\partial\hat{\varphi}^{w}}{\partial\tilde{p}}\Delta p\,dv\,.\label{eq354d} \end{equation}$$ Finally, the linearization along a concentration increment $\Delta\tilde{c}^{\iota}$ is $$\begin{equation} \begin{aligned}D\delta G\left[\Delta\tilde{c}^{\iota}\right] & =\int_{b}\delta\tilde{p}\,\frac{\partial\hat{\varphi}^{w}}{\partial\tilde{c}^{\iota}}\Delta\tilde{c}^{\iota}\,dv\\ & +\overline{\mathcal{V}}\int_{b}\delta\tilde{p}\left(1-\varphi^{s}\right)\frac{\partial\hat{\zeta}}{\partial\tilde{c}^{\iota}}\Delta\tilde{c}^{\iota}\,dv\\ & +\sum_{\gamma}\nu^{\gamma}\int_{b}\delta\tilde{c}^{\gamma}\left(1-\varphi^{s}\right)\frac{\partial\hat{\zeta}}{\partial\tilde{c}^{\iota}}\Delta\tilde{c}^{\iota}\,dv \end{aligned} \,.\label{eq354e} \end{equation}$$