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).