Consider a bivariate spatial process $\Yvec(\svec) \equiv ( Y_1(\svec), Y_2(\svec))'$, $\svec \in G$, where we refer to $G \subset \mathbb{R}^2$ as the geographic domain. The cross-covariance functions are \[ C_{ij,G}(\svec, \uvec) = \cov(Y_i(\svec), Y_j(\uvec)); \quad i,j = 1,2; \quad \svec, \uvec \in G, \] where $\{ C_{ij,G}(\cdot \, , \cdot) \}$ are nonstationary and asymmetric (cross-)covariance functions. Under warping, they can expressed as simpler, usually stationary and symmetric (cross-)covariance functions, $\{ C_{ij,D}(\cdot \, , \cdot) \}$, on a deformed space $D$.

First, consider two warping functions $\fvec_1(\cdot), \fvec_2(\cdot)$ such that $\fvec_1, \fvec_2: G \to D$. In this case, \begin{equation*} C_{ij, G}(\svec, \uvec) = C_{ij, D}(\fvec_i(\svec), \fvec_j(\uvec)) = C^{o}_{ij, D}({\fvec_i(\svec) - \fvec_j(\uvec)}); \quad i,j = 1,2; \quad \svec,\uvec \in G, \end{equation*} where $\{ C^{o}_{ij,D}(\cdot) \}$ are stationary and symmetric (cross-)covariance functions on the warped domain $D$.

In practice, any asymmetry present is likely to be simple and dominated by global shifts and rotations. To model asymmetry, we thus propose expressing each $\fvec_i(\cdot), i = 1,2$, as a composition of a shared warping function $\fvec(\cdot)$, and an aligning function $\gvec_i(\cdot)$. That is, for $i,j = 1,2$, we let \begin{align*} C_{ij, G}(\svec, \uvec) &= C_{ij, D}(\fvec \circ\gvec_i(\svec), \fvec \circ \gvec_j(\uvec)) \nonumber \\ &= C^{o}_{ij, D}(\fvec \circ \gvec_i(\svec) - \fvec \circ\gvec_j(\uvec)); \quad \svec,\uvec \in G. \end{align*} We can choose the aligning functions $\gvec_i(\cdot), i = 1,2$, to be simple transformations that are commonly used to align spatial fields and which can include translations and rotations [@wiens2019surface], and the shared warping functions $\fvec(\cdot)$ to be composed of deep injective warping functions, which are able to introduce complex nonstationary behaviour (Zammit-Mangion et al. 2019).