We adapt the viscosity method introduced by Riviere (Publ Math Inst Hautes Etudes Sci 126:177-246, 2017. https://doi.org/10.1007/s10240-017-0094-z) to the free boundary case. Namely, given a compact oriented surface Sigma, possibly with boundary, a closed ambient Riemannian manifold (M-m, g) and a closed embedded submanifold N-n subset of M, we study the asymptotic behavior of (almost) critical maps Phi for the functional E-sigma (Phi) := area(Phi) + sigma length(Phi|(partial derivative Sigma)) + sigma(4) integral(Sigma) |parallel to(Phi)|(4) vol(Phi) on immersions Phi : Sigma -> M with the constraint Phi(partial derivative Sigma) subset of N, as sigma -> 0, assuming an upper bound for the area and a suitable entropy condition. As a consequence, given any collection F of compact subsets of the space of smooth immersions (Sigma, partial derivative Sigma) -> (M, N), assuming F to be stable under isotopies of this space, we show that the min-max value inf max area(Phi) A is an element of F Phi is an element of A is the sum of the areas of finitely many branched minimal immersions Phi((i)) : Sigma((i)) -> M with Phi((i)) (partial derivative Sigma((i))) subset of N and partial derivative(nu)Phi((i)) perpendicular to TN along partial derivative Sigma((i)), whose (connected) domains Sigma((i)) can be different from Sigma but cannot have a more complicated topology. Contrary to other min-max frameworks, the present one applies in an arbitrary codimension. We adopt a point of view which exploits extensively the diffeomorphism invariance of E-sigma and, along the way, we simplify and streamline several arguments from the initial work (Riviere 2017). Some parts generalize to closed higher-dimensional domains, for which we get an integral stationary varifold in the limit.