Markov modulation is versatile in generalization for making a simple stochastic model which is often analytically tractable to be more flexible in application. In this spirit, we modulate a two dimensional reflecting skip-free random walk in such a way that its state transitions in the boundary faces and interior of a nonnegative integer quadrant are controlled by Markov chains. This Markov modulated model is referred to as a 2d-QBD process according to Ozaw (2013). We are interested in the tail asymptotics of its stationary distribution, which has been well studied when there is no Markov modulation. Ozawa studied this tail asymptotics problem, but his answer is not analytically tractable. We think this is because Markov modulation is so free to change a model even if the state space for Markov modulation is finite. Thus, some structure, say, extra conditions, would be needed to make the Markov modulation analytically tractable while minimizing its limitation in application. The aim of this paper is to investigate such structure for the tail asymptotic problem. For this, we study the existence of a right subinvariant positive vector, called a superharmonic vector, of a nonnegative matrix with QBD block structure, where each block matrix is finite dimensional. We characterize this existence under a certain extra assumption. We apply this characterization to the 2d-QBD process, and derive the tail decay rates of its marginal stationary distribution in an arbitrary direction. This solves the tail decay rate problem for a two node generalized Jackson network, which has been open for many years.