If I understand correctly, you have x_i as input and two different GP evaluate on part of x_i which gives z_i and y_j right? You can partition x_i into two part and model two GP within a pm.Model() block. Moreover, you can specify \theta one and input them into each GP if they are the same for both GP.
Not sure I understand why you said \theta is not straight forward to model - are they regular covariance function parameters like lenght scale etc?