There are lots of function optimising routines that could be applied, based on the description so far. Random search, grid search, hill-climbing, gradient descent, genetic algorithms, simulated annealing, particle swarm optimisation are all possible contenders that I have heard of, and I am probably missing a few.
The trouble is, starting with next to zero knowledge of the black box, it is almost impossible to guess a good candidate from these search options. All of them have strengths and weaknesses. To start with, you seem to have no indication of scale - should you be trying input parameters in any particular ranges? So you might want to try very crude searches through a range of magnitudes (positive and negative values) to find the area worth searching. Such a grid search is expensive - if you have $k$ dimensions and want to search $n$ different magnitudes, then you need to call your black box $n^k$ times.
This can be done in parallel though, and given you are confident that the function is roughly unimodal, you can start with a relatively low number of n (maybe check -10, -1, 0, +1, +10 for 15625 calls to your function taking roughly 8 hours 40 mins using 5 boxes). You may need to repeat with other params once you know whether you have found a bounding box for the mode or need to try yet more values, so this process could take a while longer - potentially days if the optimal value for param 6 is more like 20,000. You could also refine more closely, once you have a potential mode you might want to define another grid of values to search based around. This basic grid search might be my first point of attack on a black box system where I had no clue about parameter meaning, but some confidence that the black box output had a rough unimodal form.
Given the speed of response you should be storing all input and output values in a database for faster lookup and better model building later. No point repeating a call taking 10 seconds when a cache could look it up in 1 millisecond.
Once you have some range of values you think that a mode might be in, then it is time to pick a suitable optimiser.
Given the information so far, I would be tempted to run either more grid search (with a separate linear scaling between values of each param) and/or a random search, constrained roughly to the boxes defined by the set of $2^6$ corner points around the best result found in initial order-of-magnitude search.
At that point you could also consider graphing the data, to see if there is any intuition about which other algorithms could perform well.
With the possibility of parallel calls, then gradient descent might be a reasonable guess, because you can get approximate gradients by adding a small offset to each param and dividing difference that causes in the output by it. In addition, gradient descent (or simple hill climbing) has some chance of optimising with less calls to evaluate the function than approaches that rely on many iterations (simulated annealing) or lots of work in parallel (particle swarm or genetic algorithms). Gradient descent optimisers as used in neural networks, with additions like Nesterov momentum or RMSProp, can cope with changes in function output "feature scale" such as different sizes and heights of peaks, ridges, saddle points.
However, gradient descent and hill climbing algorithms are not robust against all function shapes. A graph or several of what your explorations are seeing may help you to decide on a different approach. So keep all the data and graph it in case you can get clues.
Finally, don't rule out random brute-force search, and being able to just accept "best so far" under time constraints. With low knowledge of the internals of the black box, it is a reasonable strategy.
