Fitting various environmental variables to an NMDS

Question

I have two sets of environmental variables (e.g. river flow and river temperature statistics). I would like to assess which of these explains the highest amount of variance of a community community (e.g. response of fish communities) that has been modeled through an NMDS. I.e. from the NMDS ordination, I would like to fit two sets environmental variables to the ordination space. Using this example, it should assess if flow or temperature exert the greatest control over the fish community. I have seen different ways of doing this, such as envfit and ordisurf functions in the vegan package of R or even with the sort of setup you get with an RDA and CCA (albeit with environmental vectors fitted to the ordination rather than having the axes constrained by each independent variables - as is needed with an NMDS).

^{- Note: NMDS = non-metric multidimensional scaling; RDA = redundancy analysis; CCA =canonical correspondence analysis.}

With my responses being non-linear, is there a function that can do this in a non-linear format that accounts for various explanatory variables? I have seen a question elsewhere NMDS and variance explained by vector fitting that suggests ordisurf (vegan package), but I would like to apply a nonlinear fitting function to multiple explanatory variables if possible, rather than the single variable used in ordisurf.