Evaluating uncertainty in solute transport under nonstationary flow conditions is a computationally demanding task. This is particularly true for cases with a two-point covariance function of log conductivity depending on the actual positions of the points rather than their distance vector. These may occur when the geological formation exhibits a trend. Nonstationarity can also be the result of uncertainty in the trend parameters of the mean log conductivity value, or it may originate from conditioning of the log conductivity field to measurements of, for example, head or conductivity. We present an efficient numerical method for evaluating the variance of travel time in such formations. We cover cases in which the nonstationary covariance functions are constructed from stationary counterparts, either by scaling functions or by summation with nonstationary functions resulting from marginalization or conditioning. We apply a matrix-based first-order second-moment (FOSM) method for uncertainty propagation, using the continuous adjoint-state method for coupled systems to evaluate the sensitivity matrix. The resulting matrix-matrix multiplications are accelerated by fast Fourier transformation (FFT) techniques after periodic embedding of the covariance matrices referring to the stationary counterparts. The combination of these methods makes it possible to compute the travel time uncertainty in domains discretized by several hundred thousand log conductivity values on standard personal computers within a reasonable time-frame. For demonstration, we apply the method to a binary medium and a medium exhibiting a continuous trend in the covariance function. In the latter application we also demonstrate the effects of marginalization and conditioning.