For multivariate functional data recorded from a sample of subjects on a common domain, one is often interested in the covariance between pairs of the component functions, extending the notion of a covariance matrix for multivariate data to the functional case. A straightforward approach is to integrate the pointwise covariance matrices over the functional time domain. We generalize this approach by defining the Fréchet integral, which depends on the metric chosen for the space of covariance matrices, and demonstrate that ordinary integration is a special case where the Frobenius metric is used. As the space of covariance matrices is nonlinear, we propose a class of power metrics as alternatives to the Frobenius metric. For any such power metric, the calculation of Fréchet integrals is equivalent to transforming the covariance matrices with the chosen power, applying the classical Riemann integral to the transformed matrices, and finally using the inverse transformation to return to the original scale. We also propose data-adaptive metric selection with respect to a user-specified target criterion, such as fastest decline of the eigenvalues, establish consistency of the proposed procedures, and demonstrate their effectiveness in a simulation. The proposed functional covariance approach through Fréchet integration is illustrated by a comparison of connectivity between brain voxels for normal subjects and Alzheimer's patients based on fMRI data.